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Fluctuating Lift and Drag Acting on a 
Cylinder in a Flow at Supercritical 


Reynolds Numbers' 


. &. 


FUNG* 


California Institute of Technology 


Summary 


The fluctuating lift and drag acting on a circular cylinder in a 
flow of an incompressible fluid at large Reynolds numbers were 
measured. Data on the root-mean-square values of the lift and 
drag coefficients, the extreme values of these coefficients, and 
their power spectra at various Reynolds numbers are presented. 


Symbols 


projected area of cylinder = length X diameter, ft.* 

lift/(qgA ) = lift coefficient 

drag/(gA) = drag coefficient 

drag force, a force component parallel to the velocity 
vector of the undisturbed flow, Ibs. 

cylinder diameter, ft. 

normalized power spectrum, see Eqs. (3), (4) 

lift force, a force component normal to the plane of 
the cylinder axis and the velocity vector of the 
undisturbed flow, Ibs. 

scale of turbulence, ft. 

mass, slugs 

frequency, cps 

power spectrum 

pU?/2 = dynamic pressure of flow, lb. /ft.* 

Reynolds number, based on cylinder diameter 

nd/U = Strouhal number 

time 

velocity of undisturbed flow, ft./sec. 

4 elastic response, see Eqs. (5), (13) 
Z(iw) = mechanical impedance, see Eqs. (5), (14) 


Received September 14, 1959; and also presented at the 
Hypersonic Flow Theory Session, IAS 28th Annual Meeting, 
N.Y., Jan. 25-27, 1960. 

+ This paper is based on two earlier technical reports of the 
Space Technology Laboratories, references 1 and 2. The study 
was instigated and supervised by Dr. Millard V. Barton, Head of 
the Engineering Mechanics Department, STL, with major con- 
tributions made by the following Members of the Technical Staff 
of the STL: Messrs. Robert Hutton, Joseph Hirsch, Robert 
Ryder, Harold Kirsch, and Edward Wenzinger. 

* Professor of Aeronautics; also Consultant, Space Technology 


Laboratories, Inc. 


= ratio of damping to the critical damping 
density of fluid 
viscosity ot fluid 
kinematic viscosity of the fluid 
time lag 
a reference spectrum, Eq. (10 
a correlation function, Eq. (6 
frequency, rad./sec 
= a reduced frequency, Eq. (9) 


(1) Introduction 


| a eanionnqant ABOUT the fluctuating lift and drag act- 
ing on a cylinder is of considerable practical interest 


in aeroelasticity, as well as in a basic understanding of 
fluid mechanics. A great deal of discussion about the 


wind-induced oscillations of smokestacks, missiles, etc., 
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Fic. 1. Mean drag coefficient and Strouhal number versus 
Reynolds number for a circular cylinder. Data Sources: Relf, 
R&M 917. Roshko, NACA TN 2913. Kovasznay, Proc. Roy 
Soc. 1949. Stack, NACA ACR. Delany and Sorenson, NACA 
TN 3038. Gé6ttingen, Ergebnisse AVA 1923. Drescher, Z. f 
Flugwiss. 4, 1956 
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Rear view of model looking toward entrance section 
of wind tunnel. 





Fic. 2(b). Model with magnesium sleeve removed. 


is found in the literature, but the basic aerodynamic 
data are lacking, particularly at higher Reynolds num- 
bers which are of practical interest for the larger 
installations. It is the purpose of the present article 
to furnish some information about circular cylinders in 
flows at supercritical Reynolds numbers. The data pre- 
sented here, concerning both a stationary cylinder and 
a cylinder subjected to forced oscillations of specific 
amplitudes and frequencies, contribute to the under- 
standing of the vibratory behavior of large missiles 
or stacks in wind. However, it was also found that the 
axial correlation of the fluctuating forces are important, 
and that significant changes can be induced by end 
effects and by small axial geometrical perturbations. 


(1.1) The Critical Reynolds Number 


For a viscous, incompressible fluid, every dimension- 
less parameter characterizing the flow around geo- 
metrically similar bodies is a function of the Reynolds 
number :* 


Udp/yp = Ud/v (1) 


* A crude approximation is R = (U in mph) (d in ft.) & 10+. 
The similarity conditions must include the cylinder-tunnel con- 
figuration (blockage and end effects), surface roughness, free- 
stream turbulence level and turbulence structure, nonuniformity 
in the flow (as in natural wind) or the presence of free surface (as 


R= 


in water). 
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where U is the speed of flow, d is a characteristic length, 
pis the density of the fluid, u is the viscosity, and v=/p 
is the kinematic viscosity of the fluid. For a circular 
cylinder the Reynolds number will be based on the 
cylinder diameter. Fig. 1 shows the variation of the 
mean drag coefficient and the Strouhal number with the 
Reynolds number. It is seen that the drag coefficient 
suffers a sudden drop when & is somewhere between 
100,000 to 500,000; the exact value of this ‘‘transition”’ 
Reynolds number zone, or the region of “‘critical’’ 
Reynolds numbers, &,,, depends on the surface rough- 
ness, wind tunnel turbulence structure and level, etc. 
The separation point of the boundary layer moves 
rearward on the cylinder in this critical Reynolds num- 
ber region. The wake near the cylinder for a flow with 
R< R., has a clear periodic structure, with a dominant 
frequency (at each side of the cylinder) expressed non- 
dimensionally as the Strouhal number S. The wake for 
a flow with R > R., is much more turbulent, and it is 
necessary to speak of a power spectrum rather than a 
single dominant frequency. Henceforth we _ shall 
speak of the Reynolds number regions less or greater 
than the transition zone as “subcritical” or “‘super- 
critical,’ respectively. 


(1.2) Review of Previous Findings 

The literature on the subject of vortex shedding in 
a flow around a cylinder has been reviewed in references 
3, 4, and 5, where extensive lists of references can be 
found. The majority of the literature concerns the 
wake development at subcritical Reynolds numbers. 
Only the works of Relf and Simmons (1924),® Delany 
and Sorensen (1953),7 and Drescher” give experimental 
data on the shedding frequency in the supercritical 
Reynolds nuinber ranges. Other than the mean drag 
component, the force component that act on the 
cylinder in supercritical flow have never been measured 
directly. It is well known, however, that there act 
on the cylinder an oscillatory lift force of the same 
order of magnitude as the mean drag and an oscillatory 
drag force superposed on the steady mean value. 
In 1879, one year after publication of Strouhal’s basic 
paper,’ Lord Rayleigh® published the observation on the 
vibration of a string under the stimulus of a chimney 
draught that, ‘“‘contrary to the opinion generally ex- 
pressed, the vibrations are effected in a plane perpendic- 
ular to the direction of the wind.’ In recent years, 
such transverse oscillations have been found to occur 
on smokestacks, submarine periscopes, oil pipe lines, 
antennas, etc., and have attracted much attention from 
engineers. (See reference 5, Chapters 2 and 9; also, 
references 10-14.) 

The available information on the fluctuating lift will 
now be reviewed. From a study of the pressure 
distribution on a circular cylinder after an impulsive 
motion, starting from rest and increasing to a constant 
velocity, Schwabe" found that when the Reynolds 


{1 The pressure was found from the velocity and the velocity 
was found from photographs of the flow. 
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gth, . number (based on the cylinder diameter) reached 735, a 
M/ p periodic lift coefficient of 0.45 and a drag coefficient of 
ular | 1.09 were attained. Bingham, Weimer, and Griffith’ 
the | determined the density and pressure distribution in a 


the — flow around a cylinder in a shock tube. At a Reynolds 
the number of about 76,000, when the first incident shock 
ient had passed ahead of the cylinder by a distance of 
yeen about 20 cylinder diameters, a nearly steady wake 
on” pattern developed at which an oscillatory lift coefficient 


CONNECTED TO 
“MG SLEEVE 





J §GI0 DISK 


“ 060 (REF) 




















eal” with amplitude of order 0.9 was obtained. 
igh- | A different measurement was made by McGregor," 
etc. | who used a microphone inside a cylinder as a pressure 
yves | transducer. Since the pressure was measured at 
um- various angular positions sequentially, not simul- S 
vith taneously, it was necessary to make certain assump- Ss 
ant | tions with regard to the phase relationship of the pres- é | 
on- | sure fluctuations at different points on the cylinders in ’ 
‘for } order to derive the resultant lift and drag forces. Mc- 
it is Gregor estimated that the amplitude of the lift co- 
na | efficient has nearly constant value of 0.6 in a Reynolds 
hall number range of 45,000 to 130,000 (subcritical). In 
ater the same Reynolds number range the oscillatory com- 
per- ponent of the drag coefficient was estimated to decrease 
| from about 0.076 to 0.044. Recently, Phillips'’ com- 
puted, from Kovasznay’s’® experimental data on the ‘ 
velocity distribution in the wake, that at Reynolds S 
rin numbers between 40 and 160 the amplitude of the lift y 
ces coefficient is approximately 0.76, and that of the g 3s 
be fluctuating drag is about 0.076. =| 5 
the The most direct measurement of the fluctuating forces c 
ers. was made by Drescher® in water tunnels in the ¢ 
any Reynolds number range 10‘ to 1.13 X 10°. The tran- “ 
ital sient pressure distribution was recorded by multiple - 
ical manometers in 12 orifices equally spaced around a - 
rag | circumference at the midsections of the models. It rn = ¢ 
the was found that the amplitude of the oscillatory lift co- [ ~y i 
red efficient varies with time and lies between 0.6 to 1.3; ' 
act similarly, the steady component of the drag coefficient 
me lies between 0.95 to 1.2; and the amplitude of the oscilla- 
ory tory component of the drag coefficient lies between 0.06 | 
ue. to 0.3. It was shown that at higher Reynolds numbers, bd 
SiC approaching the transitional regime, the fluctuating ° 
the forces are nonharmonic; they appear to come in suc- 
1ey cessive bursts. 
ex- This exhausts the published information on the 
lic- fluctuating lift and drag coefficient. 
rs, 
~ur (1.3) Outline of the Present Experiment 
es, The purpose of the present investigation is to meas- 
om ure the lift and drag forces acting on a cylinder 
so, in a flow at supercritical Reynolds numbers (of 
order 0.3 X 108 to 1.4 X 108). The cylinder model 
vill was so rigid that when it was tested in the fixed-end 
ire configuration the elastic deformation was quite negligible 
ive (the maximum static deflection under the most severe SS nn ———_ 
int condition was of order 0.005 in., the diameter being 
lds 12.65 in.). The ends of the cylinder, however, could 
also be connected to a “‘shaker’’ so that a parallel 
ity motion across the wind with variable amplitude and 
frequency could be imposed. With the inertia forces 
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Fic. 4. The shaker. 


induced by the motion of the cylinder cancelled by a 
set of mass-balances, the aerodynamic lift and drag 
could be determined. 

Since in the supercritical Reynolds number range the 
lift and drag forces are nonharmonic, it was necessary to 
determine the statistical characteristics of these forces. 
In the present experiments the lift and drag force-time 
histories were recorded on magnetic tapes, from which 
spectral analyses were made. The root-mean-square 
values of the fluctuating lift and drag were determined 
directly during testing. 


(2) Experimental Arrangement 


(2.1) The Model 


The tests were made in the the 10-ft. tunnel of the 
Guggenheim Aeronautical Laboratory, California In- 
stitute of Technology. The test model had a diameter 
of 12.65 in. and a length of 6 ft. It was supported be- 
tween two vertical channels that spanned the tunnel and 
were tied to an unusually heavy floor beam beneath the 
working section. Fig. 2a shows the rear view of the 
model (looking toward the entrance section of the 
wind tunnel). The model center section, 22 in. long, 
was carefully machined from a magnesium casting. 
The light weight of this section (effective weight, 7 lbs.) 
enabled us to design strain-gage transducers of high 
rigidity to ensure a sufficiently high resonance fre- 
quency for the force-measuring system (the lower end 
of the resonance frequency spectrum of the model in the 
simply supported condition being 208, 223, 720, 1560, 
2460 cps). The aerodynamic forces acting on this mag- 
nesium sleeve were measured. 

Fig. 2b shows the model with the magnesium sleeve 
removed. The foot of a connector that connected the 
magnesium sleeve to an inner disc and one of the three 
roll-restraint flexures can be seen. The supporting 
cylinder was made of an aluminum pipe. A disc was 
connected rigidly to the outer, magnesium sleeve. One 
end of each strain-gage element was connected to the 
disc, and the other end was connected to the supporting 
cylinder (see Fig. 3). The flexures on top of the strain- 
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gage elements were used to support beams on which 
the balancing weights were attached. 

The function of the strain-gage elements can be 
discerned more clearly from the drawing in Fig. 3. 
The gages read drag andliftseparately. The magnesium 
sleeve and the strain-gage elements were connected 
with two balancing masses, one in the lift plane and 
one in the drag plane. These masses were positioned by 
screwing them on their supporting shafts in such a 
way that the fundamental component in the strain- 
gage output (of the forcing frequency) was minimized 
when the model was forced to execute sinusoidal oscilla- 
tions of various frequencies and amplitudes in still air. 
Thus the inertia forces caused by the forced oscillations 
were essentially eliminated from the strain-gage read- 
ings which, therefore, measured directly the aero- 
dynamic forces. 

The ‘‘shaker,”’ 
oscillations on the cylinder in the lift plane. This 
mechanism consisted of a variable-speed motor which 
drove a shaft below the tunnel by means of pulleys and 
belts. Two eccentric couplings, connected at each end 
of the shaft, imposed an up-and-down reciprocating 
motion in the lift plane on the model in the wind tunnel 
above and on a set of balancing weights below. The ec- 
centricity of the couplings could be set for a desired 
amplitude of motion. For stationary cylinder tests, 


shown in Fig. 4, imposed forced 


separate end supports were attached to lock the model 
positively at its ends. 


(2.2) Instrumentation 

The strain-gage transducers were specially made for 
the present test. 

The lift and drag forces and accelerations were re- 
corded on an Ampex FR-100 modular magnetic tape 
recorder with an FM system. Fig. 5 shows a block 
diagram of the instrumentation set-up. The _ root- 
mean-square values of the strain-gage signals were read 
from a Ballantine ‘““True RMS” voltmeter. In order to 
eliminate contributions to the root-mean-square values 
due to resonance frequencies of the model, a filter was 
used. 

Static calibrations were made with dead weights, and 
were recorded on the magnetic tape for every test run. 
Linearity and symmetry of the strain-gage outputs 
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Fic. 5. Instrumentation set-up. 
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Fic. 9. The steady and oscillatory components of the drag co- 
efficient for a fixed cylinder. 
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were assured to within 1 per cent in the full load range. 
Alignment in the wind tunnel was made by minimizing 
the ‘‘cross-talk’’ between lift and drag gages with respect 
to deadweight loadings. 

A dynamic calibration for the root-mean-square volt- 
meter reading was made by strapping a small eccentric- 
weight shaker on the magnesium sleeve and noting the 
corresponding voltmeter readings and the shaking fre- 
quencies. Two identical eccentric weights, symmetri- 
cally disposed and rotating in opposite directions, im- 
posed unidirectional harmonic forces on the cylinder. 
The voltmeter reading could then be converted into 
the root mean-square values of the exciting force 
These calibrations were further checked by the cali- 
brated oscillograms. 

Since the wave form was nonharmonic at large 
Reynolds numbers, a question arose as to the reliability 
of the voltmeter in giving the root-mean-square values 
No detailed investigation was made in the present test, 
but a comparison of the simultaneous readings from a 
Hewlett-Packard peak-reading voltmeter showed that 
the “‘true r.m.s.’”’ readings may differ from the peak- 
reading meter data by as much as 50) per cent, with the 
Ballantine giving higher values generally. The ratio 
(Ballantine reading)/(Hewlett-Packard reading), for 
an overwhelming majority of cases, fell between 1.1 and 
1.3. The Ballantine data were considered as more 
applicable to the present purpcse and are presented be- 
low. 


(2.3) Data Reduction 

Oscillograms were obtained from the magnetic tapes 
and recorded on Hieland recorder. Samples of these 
records are given in Fig. 6. The +7.5 per cent fre- 
quency-deviation lines were printed on the oscillogram. 
The neutral position (with no aerodynamic load) is mid- 
way between these deviation lines. A displacement of 
the signal from the neutral position indicates a load on 
the magnesium cylinder. The load-displacement re- 
lationship was linear. 

The spectral analysis was made with an automatic 
wave analyzer of the heterodyne type, Model 9020, 
Davies Laboratories, Honeywell. In general, a 2-sec. 
loop (120 in.) was made from the magnetic tape record- 
ing, and the spectrum of this loop was analyzed. In 
order to compare sampling errors, several 10-sec. loops 
were also made. Fig. 7 illustrates the amplitude versus 
frequency spectra plotted automatically by the ana- 
lyzer. In these plots a fixed (3 decibel) bandwidth 
setting of 6 cps was used. 


(2.4) Accuracy 

Since the frequency range of interest was below 200) 
cps, the static and dynamic calibration procedures 
described above were satisfactory. The probable error 
in the scale for the oscillograms was within | per cent. 
The calibrated scale for the root-mean-square voltmeter 
reading was within 2 per cent. The effect of the non- 
harmonic waveform was unknown. 
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Fic. 10(b). 








The major contribution to the probable error in the 
final results was due to the randomness of the phe- 
nomenon being studied. The needle of the voltmeter 
fluctuated up and down, even after the time constant 
was lengthened. Correspondingly, the power spectra 
varied from sample to sample. 
study of a Gaussian noise, Magness”! has shown that 
the standard deviation of a measured power spectrum 
is of order 14/B7, where B is the effective bandwidth 
of the analyzer and 7 is the length of the sample-tape 


From a theoretical 


loop. As shown in a paper by Barr and Morrison,” a 
setting of 6 cps and a 2-sec. loop corresponds to a varia- 





(Len?) 2 


[ 2 
| Reference Curve 6(n) -< 2(143n*) 
| 
| 


Test Point 305 





| | 
! 

















q = 41b/fr? 
U = 60.4 ft/sec. 
a i ee \ | R = 0.33x 10° 
&, Ps \ | | } 
‘y coe ee | | T ] ] 1 
v \ | | | 
‘ 
& a 
| | | 
. P| 
a, 4+ SS eS aE 
. " T T yr re 4 
V j 
§ } 
g N\ | 
e | \\ | 
$3 . a 4 — Se — EE —EEEE 
‘ \ T T° 
H \ | 
: \ 
\ | | | 
Hi \ | | 
"2 + * T — + + = Ge | 
| \ | | 
7 \ | | | 
. \ | 
3 \ | 














Nondimensional Frequency, S = nd 
U 


Fic. 10(a). Normalized power spectrum for the lift force at 


Reynolds number 330,000. 














‘eee ne | 
| | a 
| Reference Curve (nm) = e ft 
| | | 4 en’) 
6 | | | Test Point 306 
| | q = 5 lb/ft? 
| U = 67,5 ft/sec 
se] A } R = 0,38 x 10° 
i ’\ 
| on —— 
| | | | | 
: + + | — | 
| 









L || 








Normalised Power Spectrum of Lift Force, 











0,3 


4 
Nondimensional Frequency, S = “s 


0.4 0.5 ° 0.7 


Normalized power spectrum for the lift force at 
Reynolds number 380,000. 


FLUCTUATING LIFT AND 





Fic. 10(d). 








DRAG 








However, since the statistical dis- 


tion of 16 per cent. 
tribution of the fluctuating forces was unknown, no re- 
liable estimate of the sampling error could be made. 
Further discussion will be given to this later on. 

An attempt at an estimation of the variation in the 
root-mean-square values of lift and drag was made by 
noting the range of fluctuations of the voltmeter needle 
during the test. In general, the variation from the 
mean did not exceed +7 per cent, but on rare occasions 
it was as high as + 15 per cent. 

For the forced oscillation tests there were two addi- 


tional sources of errors. First, there might be errors 
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in the shaker amplitude settings, which were made by 
matching a pair of cross-serrated couplings. The pitch 
of the serration was 1/8-in. Any mismatch in the 
serrations would cause an error in the amplitude and 
any difference in amplitudes at the two ends of the 
model would cause an undesirable rocking of the model. 
Assuming that in the worst case a mismatch amounted 
to 1/4 of the pitch of the serration and that the two 
ends of the model mismatched in opposite directions, 
then a rocking motion with an amplitude of about 3.5 
min. was possible. In general the mismatch was con- 
siderably less than this, but the effect was noticeable. 
The second source of disturbances came from the vibra- 
tions of the shaker linkage and the vertical channels 


which supported the model. These vibrations were 


= 30 lb/ft” 


i | - : | | > 
| ‘i 2 
| | | Reference Curve (nm) = eee) 
| d 2,2 
| | (l+n") 
} 
‘+--+ | = Test Point 311 
| | | > 
| 


| q 














a 
‘ U = 165 ft/sec 
3 R = 1,00 x 10° 
* 
5 +- = | TN r _ 
lee | 
s 
3 
% 
54 —_+—— + + + “4 
2 
Fa 
7) 
& | | 
a | 
be | | 
$3 + —- + } 4 
° | | 
= 
v | | 
* 
F } 
2 = 4 + } 4 1 
& \ 
Zz ‘ 














» d 
Nondimensional Frequency, S = 7 


Normalized power spectrum for the lift force at 
Reynolds number 1,000,000. 


Fic. 10(e) 


| 
| 
| 
| 
















Reference Curve $(n) = 


Test Point 31? 
q = 40 Ibe/ft’ 
U = 191.0 ft/sec 
R = 1,13 x 10° 


~ te 





| 

| 

| 

| | | 
; ne ieee 
| | 

| | 























Nondimensional Frequency, S = — 
U 


Normalized power spectrum for the lift force at 
Reynolds number 1,130,000. 


Fic. 10(f). 


SCIENCES 1960 


-NOVEMBER, 


excited by the shaker. The frequencies of the bending 
vibrations of the vertical channels in the drag direction 
were too low to permit reliable measurements of the 
oscillatory component of the drag force in the forced 
oscillation tests. 


(3) Results and Discussion 


(3.1) The Lift Coefficient C, on a Stationary Cylinder 


The lift force, transverse to the direction of flow, 


has a steady mean value of zero. Fig. 8 shows the 
amplitude of the lift coefficient. The root-mean-square 
value remains nearly constant for large Reynolds num- 
bers, but increases as the Reynolds number decreases into 
critical Reynolds 


the transition region. The exact 
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number for the test model was not determined, but it 
appeared to be somewhat lower than those obtained by 
other authors on smaller models. The peak values of 
C, shown in this figure (the largest value on the record) 
were obtained from the oscillograms (recorded over a 
2-sec. period). 

It should be recalled that for subcritical Reynolds 
numbers an amplitude of C, ranging from 0.45 to 0.9 
has been given by various authors. 


(3.2) The Drag Coefficient C, on a Stationary Cylinder 

The root-mean-square values and the peak values of 
the oscillatory component of the drag coefficient were 
obtained in the same manner as the lift coefficient (see 
Fig. 9). The scatter is larger than that of the C,, prob- 
ably in part because of a larger contribution by a reso- 
nance condition around 200 cps. The mean value, or 
the steady component, was obtained from the oscillo- 
gram and has a value of approximately 0.25 for super- 
critical Reynolds numbers. This value agrees reason- 
ably well with Delany and Sorensen’ at R=0.5 X 105, 
but is considerably smaller than the Cp given in refer- 
ence 7 at higher Reynolds numbers. For example, 
reference 7 gives a range of Cp varying from 0.29 to 
0.45 at R=1.4 X 108 Delany and Sorensen dis- 
cussed the effect of leakage at the ends of their model. 
In our test, no attempt was made to seal the ends of the 
model. The wide variation of Cp in reference 7, as 
well as the difference between their result and ours, 
can probably be explained by differences in end condi- 
tions and the free stream turbulences, but this is not 
yet verified. 
(3.3) The Power Spectrum of the Lift Force on a Stationary 
Cylinder 

Power spectra for the lift force at Reynolds numbers 
between 0.335 X 10° and 1.39 X 10® are plotted in 
normalized form in Fig. 10. The frequencies are ex- 
pressed in terms of the nondimensional frequency 5, 
the “Strouhal number’’: 


S = nd/U (2) 
where m is the frequency in cycles per second, d is the 
cylinder diameter in feet, and U is the speed of un- 
disturbed flow in feet per second. The power spectrum 
F(S) for the lift force is so normalized that* 


f F(S)\dS = 1 (3) 


The mean square lift, on assuming no phase difference 
in the axial direction, is therefore given by 


L? = @AyCe f F(S)dS (4) 


where g is the dynamic pressure of the flow, A is the 
projected area of the cylinder, and C,;? is the mean 
square value of the lift coefficient. 


* The spectra in Fig. 10, however, are normalized with an up- 
per limit of 480 cps. 


















































7 . 4 
id @ = 40 tbe/h* 
U «= 191 &/eec 
Vv © 
a> v R = 1,128 = 10 
6 
a 
2 qeee 
u B hy © 10-sec loop 
§ te o, : 42094 D 2-sec loops 
55s x 
~ wv pe? 
< a 
4 *a5qe 
) oO 
E«4 Ls 
4 4,0 
7 td 
& vb 
= Kd 
o 
Bs $ La} 
é BSpa 
© 70 
- 4, 
a a", 
E. Soe 
- 
5 
1 4 $ 
pi, 
4 
We $ 
al AARE EE PST T 
0 0.1 0.2 0.3 04 0.5 0.6 0.7 





Nondimensional Frequency, S = ad 
U 


Fic. 11. Variations of the lift spectra with samples. A compos- 
ite plot for five 2-sec. samples and one 10-sec. sample. 


If the lift force acts on an elastic structure whose 
mechanical impedance is Z(iw). and if the induced mo- 
tion is sosmall that the aerodynamic action is unaffected, 
then the mean square value of the elastic response of 
the structure, y, can be computed from the expression 


-—? FS 
f= (gA)*C,? - dS (5 
y gA)*C1 f |Z (ica) |? C (5) 


(see p. 297 of reference 5, or reference 23). The im- 
portance of this application is a principal reason for 
measuring the power spectrum of the aerodynamic 
forces. An example is given in the Appendix. 

Since the lift force at high Reynolds number is a 
random function of time, a correlation function ¥(7) can 
be defined as the limiting value of 


¥(r) = Lit) L(t+r) (6) 


averaged over a long period of time. The power spec- 
trum and the correlation function are, when properly 
normalized, Fourier cosine transforms of each other (p. 
289 et seq., reference 5). If the form of the correlation 


function were 
Vir) = W(O)e "1 — (rU/2D)) (7) 
where / is a characteristic length called the “scale of 


turbulence,’ then the corresponding power spectrum 
would assume the form 


p(w) = p(0){(1 + 3n?)/(1 + n*)?] (8) 
9 = wl/U = 2x(l/d)S (9) 


where ¥((), p(0) are the values of ¥(7) and p(w), respec- 
tively, at r=w=0. These are the characteristic func- 
tions for the transverse velocity components in wind 
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Fic. 12(a). Normalized power spectrum for the drag at 
Reynolds number 1,390,000. 


tunnel turbulences at nearly isotropic conditions (p. 296, 
reference 5). 
In Fig. 10, a reference curve 


(5) = 2(I/d) {{1 + 3(2a1S/d)2]/[1 + (2xlS/d)2]2} (10) 


for which 


f ¢(S)dS = 1 (11) 


is sketched for comparison with the measured spectra. 
The value of the ratio //d was taken as an average that 
seems to fit most of the curves: 


I/d = 2.4 (12) 


From Fig. 10, one is tempted to conclude that the 
power spectrum of the lift force does not vary signifi- 
cantly with the Reynolds number in the supercritical 
range. 

The character of the spectra in Fig. 10 is markedly 
different from that in the subcritical Reynolds number 
range. In the latter range one would expect, from 
Roshko’s results,‘ a large peak (in fact, a delta function) 
to appear in the power spectrum. If a universal Strou- 
hal number based on the wake width, as suggested 
by Relf® and Roshko,* remains unchanged as the 
Reynolds number passes from subcritical to super- 
critical range, one would expect a delta function to 
appear at S~0.2 in the subcritical Reynolds number 
range and at S~0.4 in the supercritical Reynolds 
number range. It is not permissible to say that our 
data did not show the existence of such peaks, but the 
size of such peaks was evidently not large. 

Fig. 11 shows an example of the variation of the lift 
spectra with statistical samples. 


(3.4) The Power Spectrum of the Drag Force ona 
Stationary Cylinder 


A few power spectra of the drag force are presented in 
Fig. 12, which are normalized with an upper limit of 
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150 cps. The steady component was eliminated so 
that each spectrum was normalized for the oscillatory 
component only. Due to the smaller signal-to-noise 
ratio, and a resonant condition of vibration in the drag 
direction at 200 cps, the accuracy of the drag spectra is 
not as high as that of the lift spectra, and the numerical 
values should not be taken too seriously. 


(3.5) The Lift Force Acting on a Cylinder Subjected to 
Forced Oscillations Perpendicular to the Flow 

Some results were obtained in oscillation tests in 
which the model was forced to execute harmonic 
oscillations of specific amplitudes and frequencies per- 
pendicular to the flow. The scope of the experiments 
was limited, however, and further studies will be needed 
to delineate fully the effects of motion on the forces 
generated in a separated and turbulent flow. 

Figs. 13-17 show the amplitude of the lift coefficient. 
The root-mean-square data were reduced according to 
expressions derived in reference 2. The peak values 
were obtained from the oscillograms, and represented 
the highest amplitudes occurring in a 2-sec. recording 
period. No corrections for disturbances due to the 
model-supporting system and the shaker were made for 
the peak values. The large scatter in peak values is 
probably caused, at least in part, by such disturbances. 

The dotted lines in Figs. 13-17 represent the average 
“stationary” cylinder values to be compared with the 
oscillating cylinder data. 

The scatter of the peak values is too broad to draw 
any conclusions. The scatter of the root-mean-square 
values is reasonably small. The general appearance of 
the data in these figures shows that the root-mean- 
square values of the lift coefficient for various oscilla- 
tion amplitudes and frequencies all fall in the same band 


of scatter. Hence it is impossible to delineate the effect 
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amplitude, 1/2 in. 


of the amplitude and frequency of the forced oscillation 
onthe lift. Therefore, a tentative conclusion is reached 
that for Reynolds numbers larger than 600,000, for 
Strouhel numbers (based on the forcing frequency) less 
than 0.12, and for a ratio of oscillation amplitude-to- 
cylinder radius of less than 1:12, the lift coefficient is 
not influenced by the forced oscillation to any large ex- 
tent. 

Fig. 18 shows a composite plot of the normalized 
power spectra for the lift at various forcing frequencies. 
The general character of the spectrum remained un- 
changed in the range of forced-oscillation Strouhal 
numbers examined here (0—0.074), for an amplitude of 
1/2 in. and a Reynolds number of 1,130,000. 

Fig. 19 shows the steady component of the drag co- 
efficient for a cylinder oscillating at an amplitude of 
1/8 in. at various frequencies. 


Concluding Remarks 


The most important feature of the forces induced on a 
circular cylinder by vortex shedding in a supercritical 
Reynolds number range, is their randomness. Recog- 
nizing this fact, we see that if the response of an elastic 
cylinder to vortex shedding at supercritical Reynolds 
numbers is to be computed, generalized harmonic 
analysis is appropriate. 

The present paper presents experimental results on 
the lift and drag coefficients in flows with Reynolds 

mbers up to about 1.4 X 10°. 

As to the effect of the motion of the cylinder itself, 
in the limited range of forcing frequencies and ampli- 
tudes studied here, the root-mean-square value of the 


lift coefficient is not influenced by the forced oscillations 
to any large extent. No drastic change in the power 
spectrum was found at these supercritical Reynolds 
numbers. We must remark that things may be quite 
different in the transitional or subcritical range. For 
example, Goldman’s*6 result seems to show a large effect 
of cylinder motion at a Reynolds number on the order of 
Lb xX 10. 

The mechanism of vortex shedding still remains 
poorly understood and defies theoretical treatment, al- 
though recent works by Roshko*! on splitter plates and 
by Phillips'* on the sound field seem to open up new 
vistas of research. It is hoped that the information on 
the supercritical flow condition and the forces acting on 
the cylinder will help toward further understanding of 
the phenomenon. 

There remain a number of questions that must be 
answered by future studies. The first is the effect of 
free-stream turbulence in the wind tunnel on the forces 
acting on the cylinder; there are indications that the 
wind tunnel turbulence level and its spectrum character- 
istics may have some important effect. The second 
is the axial correlation or the three-dimensional 
In our experiments, the model 
between each end of 
These 


character of the flow. 
has a spacing of 0.020 in. 
the center section and the outer sections. 
cracks are approximately ().9 in. deep and are not sealed. 
Wake surveys using a rake of pitot and static tubes 
showed that the wake was affected to a large extent by 
the cracks and the ends of the model. In fact, any 
small circumferential disturbances, such as a thin layer 
of tape, caused large changes in the local wake width 
and shape. Although further experiments on our par- 
ticular model showed that the forces on the central 
magnesium sleeve were not affected seriously by such 
small disturbances, the phenomenon was intriguing 
enough to warrant careful studies. Experiments by 
other aircraft companies (see references 26, 27 and some 
unpublished reports of the Martin Aircraft Company, 
Denver) have shown that the end effects of the cylinder 
can be very important in determining the amplitude of 
the oscillation of cantilevered structures. 
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Fic. 18. Normalized power spectra for lift at various forcing 
frequencies. Ratio of peak-to-peak oscillation amplitude to cyl- 
inder diameter, 1:12.6; Reynolds number, 1,130,000; forced 
oscillation amplitude, 1/2 in. 
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Appendix 
A Sample Calculation of Elastic Response 


Let a uniform circular cylinder of effective mass m, 
supported on a set of springs, be exposed to a flow of 
velocity U’ perpendicular to the axis of the cylinder. 
Let the natural vibration frequency of the spring- 
mass system in still air be w radians per second. The 
equation of the transverse motion of the cylinder is 


(d°y/dt*) + 2yuo(dy/dt) + wo?y = (1/m)L(t) (13) 


where L(t) is the force acting on the cylinder, y is the 
transverse displacement (perpendicular to the direction 
of flow), and y is the ratio of the damping constant to 
the critical damping. L(t) depends on the Reynolds 
number. If the Reynolds number is subcritical, L(f) 
is essentially periodic, and the response y(t) can be 
calculated in the usual way. (See p. 70, reference 5). 
If the Reynolds number is supercritical, L(/) is random, 
and the mean square response may be calculated (see 
reference 23 or p. 297 et seq., reference 5). However, 
the experimental results described above do not permit 
extrapolation to a cylinder of arbitrary length, because 
the correlatior >of the force along the length is unknown. 
But if we assume that the end effect is negligible, the 
flow uniform and in phase along the cylinder axis, and 
that the motion of the cylinder does not affect the aero- 
dynamic force to any significant degree, then an upper 
bound to y? may be obtained. 
In this case, the structural impedance is 


Z(iw) = (a —_ w? + 12-yaw)m (14) 


Under the assumptions listed above, we have 


L? = WAyC? [ F(S)dS = 


oe > 7 o- e 
(qA)?C;? ree | F (=) dw (15) 


and 
- d 
2 = (gA)2C,?2 
J (gA)?C, onU x 
a F(wd/22U)dw ; 
f ° | » . / le (16) 
9 mM*e*|1 — (w/wo)? + 2ty(w wo) |? 


If y < 1, an integration gives the approximation: 


] 1 ro Lt @ 
ny? = = (gA)?C;? = F ( a ) (17) 


L 2rU mMwo? Y 


We must remark that the neglect of end effects 
and axial correlation is an unjustifiable oversimpli- 
It seems that much reduction in the ampli- 
tude of the wind-excited oscillation on cantilevered 
structures can be gained by proper design of the 
geometry of the ends. We should also point out that 
Scruton ef al. have shown" that a chimney stack with a 
twelve-sided external polygonal section may have 


fication. 


“galloping instability’? and possesses certain features 
that resemble a square-section stack more than a cir- 
cular one. 
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Fic. 19. The mean values (the steady component) of the 
drag coefficient for an oscillating cylinder; forced oscillation ampli- 
tude, 1/8 in. 
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Insofar as revealing the randomness of the fluctuating lift and 
drag on a cylinder in flows of supercritical Reynolds numbers, 
these papers are in agreement with the present author. The 
papers by Buell and Ezra et al. discuss the tip effect. The papers 
by Fujino, Nakagawa, et al., give measurements of forces in towing 
tank and wind tunnel, but the statistical methods were not used 
in data evaluation and the power spectra were not given. Go- 
recki’s thesis presents an interesting account of the related prob- 
lem of temperature fluctuation in flows of Mach number between 
0.35 to 0.70, and Reynolds number of order 10°. Macovsky’s 
paper discusses the three-dimensional character of the flow on an 
ostensibly two-dimensional cylinder model. This three-dimen- 
sional character is treated in much greater detail by Humphreys. 
In this respect the theoretical investigations by Rosenhead and 
others, reviewed in reference 4, are remarkable. Finally, the ef- 
fects of artificial disturbances, such as winding of wires on the 
surface of a cylinder, are dealt with by Nakagawa,et al., and 
Weaver. 

This burst of interest is very encouraging, but it is clear that 
much work remains to be done. It is the author’s belief that the 
measurement of the correlation functions of the lift and drag 
forces along the axis of a cylinder, with and without a tip, is of 
prime importance. Such a program is being pursued at the 
California Institute of Technology by the author and Louis 
Schmidt. The theoretical background of the new Caltech ex- 
perimental program is given in a paper entitled The Analysis of 
Wind-Induced Oscillations of Large and Tall Cylindrica! Structures, 
by the author, published as Technical Report EM 10-3, Space 
Technology Laboratories, June, 1960. 
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Plastic Bending of a Work-Hardening 


Circular Plate With Clamped Edge 


CHINTSUN HWANG* 


The National Cash Register Company 


Summary 


The plastic bending of a rigid, work-hardening circular plate 
with clamped edge and under distributed transverse load is 
studied in this paper. The plate material is assumed to follow 
the yield condition of von Mises and to harden isotropically. 
The incremental stress-strain law is used. The possibility of 
partial unloading of the plate material during the loading cycle 
is investigated. The problem is represented as two sets of two 
simultaneous differential equations corresponding to the loading 
and unloading portions of the plate material. They are solved 
by the NCR 304 Digital Computer. 


Introduction 


i MODERN space-vehicle structures where extra-high 
stresses are experienced in relatively short duration, 
it is not uncommon to load the structural components 
beyond the yield point. In this state of stress all the 
metals exhibit the work-hardening effect. The stress- 
strain relation of the material becomes nonlinear, which 
in turn causes a nonlinear load-deflection relation in 
the structure. Because of the mathematical difficul- 
ties involved, the work-hardening behavior of the 
structural components has been explored only spar- 
ingly (e.g., see references 1, 2, and 3). With the ad- 
vance in the digital computer techniques in recent 
years, it seems that the vast field of the work-harden- 
ing behavior of structures can now be tackled with the 
help of the new tool. This paper, which represents 
an effort in this direction, gives the results obtained 
in solving a work-hardening plate problem through the 
use of a digital computer. 

The plate material under study is rigid, work-harden- 
ing; it follows the yield condition of von Mises and 
hardens isotropically. The incremental flow rule‘ 
is used to formulate the compatibility equation which 
is solved simultaneously with the equilibrium equa- 
tion. The boundary condition at the clamped edge of 
the circular, work-hardening plate is investigated. 
Also studied is the possible partial unloading of the 
work-hardening plate at the loading cycle of the struc- 
ture. A method to deal with cases where unloading 
does take place in portions of the structure is presented. 
Numerical data are tabulated. 


Received August 31, 1959. Revised and received January 12, 
1960. 

* Head, Commercial Systems; also, Lecturer, University of 
Southern California. 


A Work-Hardening Plate 


Consider a circular plate made of rigid, work-harden- 
ing material. The radius of the plate is RX, the thick- 
ness is 2h, which is small compared to the radius. The 
radial axis is indicated by r and the time is indicated 
by ¢. The edge of the plate is clamped. Axisymme- 
trically distributed transverse loads p(r, ¢) are applied 
to the plate so that the transverse shear force is V(r, f). 
The distributed loads are increased in a stepwise fashion 
in the work-hardening stage. 
in the plate are assumed to be essentially due to plate 


The stresses generated 


bending. The membrane stresses as well as the 
transverse shear stresses are neglected. Also neglected 
is the normal stress along the direction of the plate 
axis. 

When a rigid, work-hardening plate element is bent 
along the radial and the circumferential axes, it is 
assumed that the material normal to the initially plane 
middle surface stays normal to the deflected middle 
surface. The strain ratio along the radial and the cir- 
cumferential directions, e,:e,, remains constant along 
each normal. Suppose we have a yield curve of the 
material corresponding to plane stresses. The above- 
mentioned conditions for plate bending are satisfied 
if, along a certain normal, we represent the stresses at 
one side of the plate middle surface by a point on the 
yield curve, and represent the stresses at the other side 
of the middle surface by another point, which is dia- 
metrically opposite to the first stress point on the yield 
A state of discontinuity exists along the plate 
Based on the above assumptions, it is 


curve. 
middle surface. 
possible to use the yield curve in plane stress as the 
yield curve for plate bending with appropriate modifi- 
cation of the coordinate axes. Thus, for the case of 
plane stress, von Mises’ criterion reduces to 


=2 


a1” — o102 + oa.” = & = f(a, a2) (1) 


where a; and a are the non-zero principal stresses, & is 
the equivalent stress. Similarly, for a plate made of 
isotropic work-hardening material which obeys von 
Mises’ yield condition and a specific flow rule, the yield 
and post-yield conditions in the loading stage can be 
represented by 

F(M,, M,) = M? —- M,M,+M,2 = M* (2) 
where M, and M, are the radial and circumferential 


bending moments of the plate element, J/ is a value 
equal to or greater than the plate yield moment My = 
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M 


Mo C 





-Mo 





The bending moment vs. curvature plot in simple 
bending. 


Fic. 1. 


ooh”, oo being the yield stress of the material in simple 
tension. 

For a rigid, work-hardening plate, the moment versus 
plate curvature plot in simple bending is as shown in 
Fig. 1. When the bending moment is smaller than 
Mo, the corresponding curvature change is neglected. 
When the moment is greater than Mo, the slope of the 
moment against the curvature is assumed to be con- 
stant and is denoted by C. 

In the work-hardening range, we denote the incre- 
ment in the distributed load on the circular plate by 
Ap. The equilibrium equation for the circular plate 
corresponding to the load increment is 

= = (—AM, + AM,) — Apr (3) 

dr r . 2 
To satisfy the continuity condition for the work-hard- 
ening plate, a compatibility equation is needed. This 
equation is obtained through the use of the incremental 
flow rule. The incremental flow rule, as given by Hill,‘ 


1S 


1AM l 
(M? — 4M,M, + 4M,?) ( _ 
, 


[((4M,? — 4M,M, + 
dr 


1960 
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(@ = 1, 2, 3) (4) 


where e; are the plastic strains along the principal 
directions, o;’ are the components of stress deviation, ¢ 
is the equivalent stress, H’ is the slope of the equivalent 
stress/plastic strain curve. For the special case of plane 
stress, the flow rule can be written as 


1 of (of 0o;)do, + (Of, Oo2)daz 
4H’ Oo; f 


da, = 


(0 
je 1 Of (Of/00;)do, + (Of/O02)da2 ; 
we af On } 

where f is as defined in Eq. (1). The third strain incre- 
ment de; can be obtained from the condition that the 
sum of the strain increments in the three principal 
directions vanishes. Correspondingly, for the work- 
hardening plate, we have the following flow rule: 


1 OF (OF/0M,)AM, + (0F/0M,)AM, | 


ye F 
(6) 
ag. — 1 OF QF/OM)AM, + (OF/OM,) AM, 
*  4C0M, F 


where Ak, and Ak, are the curvature increments along 
the radial and circumferential directions, C is a plate 
constant defined previously in the paper, F is as de- 
fined in Eq. (2). 

In incremental form, the plate curvature continuity 
condition is 


dAk, 1 
“ = - (Ak, — Ak,) ( 
r 


¢ 


~] 
~— 


dr 


Combining Eqs. (6), (7), and making use of Eqs. (2), 
(3), we obtain the following compatibility equation: 


M,2)AM, + 


dr 


_ dM, . dM, 
(—M,? + 4M,M, — 4M,2)AM,] — | (—4M, + 5M,) —~ + (6M, — 4M,) —*] AM, — 


1 
. (M, — M,) ( 


dr 


(M,? — 4M,M, + 4M,2)AM,] 


The application of the above compatibility equation 
(8) to the work-hardening plate problem will be ex- 
plained later. 


The Clamped-Edge Condition 


The problem of incipient plastic flow for a circular 
plate with clamped edge and under uniformly distrib- 
uted load has been solved by Hopkins and Wang? using 
the isocline method. 





r 
(2M, — M,)(dM,/dr) — (M, — 2M,)(dM, dr) 


aM, 4, oe) AM, + [(-—2M,? + 5M,M, — 2M,2)AM, + 


M? — M,M, + M,? 


l y ‘ 
= (2M,2 — 5M,M,+ 2M,2)Apr (8) 


Under von Mises’ yield condition, the load corre- 
sponding to plastic flow based on this method is pp = 
12.5 M,/R?. The same problem solved by the digital 
computer yields a load corresponding to incipient plas- 
tic flow pp = 12.551807 Mo/R*®. The plate moment 
distribution as of this state is given in Table 1 and Figs. 
2-4. In solving the problem, it is assumed that 


Poisson’s ratio u is equal to 1/2. With further increase 
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in the distributed load, the plate material enters the 
work-hardening range. 

Corresponding to the incipient plastic flow, Eason® 
established the velocity field of the circular plate under 
distributed load in association with the stress field. 
In his study on a plate with clamped edge, it was found 
that a singularity in the curvature rate exists at the 
At that location, the slope of the deflected 
As a result, the circum- 
At the same 


edge. 
middle surface is undefined. 
ferential curvature rate is also undefined. 
time, the radial curvature rate approaches infinity. 
An arbitrary value is assumed for the slope of the de- 
flection curve at the clamped edge. The deflection 
curve is represented by a group of integral terms which 
are suitable for numerical evaluation. 

For a rigid, work-hardening circular plate with 
clamped edge, no deflection of the plate is encountered 
until the state of incipient plastic flow is reached. In 
the work-hardening range, it is expected that the 
plastic flow is contained in the plate. If w(r, ¢) repre- 
sents the deflection of the plate middle-surface, we 
have the following expressions for the curvature incre- 
ments of the plate along the radial and circumferential 
directions: 


Ak, = A(d?w/dr?) l (9) 
Ak, = A[(1/r)(dw/dr) |f 


Mo 

Mo 
2.0 
1.8 


of 


t 
or Nw 


£210 





At the clamped edge, normally we expect the curva- 
ture increment Ak, to vanish because the slope (dw/dr) 
vanishes. At the same time, the curvature increment 
Ak, is a finite quantity. Referring to the flow rule (6), 
we see that the above conditions of Ak, and Ak, at the 
clamped edge are satisfied if 


OF/OM, = —M, + 2M, = 0 (10) 


Eq. (10) is represented by the straight line marked 
r/R = 1 in Fig. 2. Substituting Eq. (10) into Eq. 
(8), it is found that the derivative (dAM,/dr) ap- 
proaches infinity at the clamped edge. 

In this paper, a numerical integration method is 
used to solve Eqs. (3) and (8) simultaneously to ob- 
tain the moment distribution of the work-hardening 
plate. The integration is carried out in a stepwise 
fashion from the center of the plate toward the clamped 
Some approximation is involved in each step- 
This approximation is clearly indi- 


edge. 
wise integration. 
cated in Fig. 4 where the computed circumferential 
bending moment M, is plotted against the radial axis. 
At the clamped edge (r/R = 1), the slope of the M, 
curve should approach infinity as explained above, 
while in the actual plot this is not true. In order to 
attain a predictable limit of possible error introduced 
by the extrapolation process at the clamped edge, 
several runs are made for each integration process from 


Mr 


"1.6-14 -12-10 -8-6-4-2 0 2 4 6 8 101214 16 1820 M, 


Fic. 2. 


The radial and circumferential bending moment distribution of a clamped circular plate of rigid, work-hardening material 


under uniformly distributed transverse load. 
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Fic. 8. The radial bending moment distribution of the rigid, 
work-hardening plate. 


the center of the plate toward the clamped edge. The 
detail of the process will be explained in the next section. 


Unloading of the Plate Elements 


In some cases, while there is an increase in the uni- 
form load Ap, it is possible that unloading occurs in 
part of the plate material. For such portion of the 
plate, the equilibrium equation (3) still holds, whereas 
the compatibility equation (8) does not apply. In this 
case, even though the curvature change corresponding 
to unloading approaches zero and can be ignored (see 
Fig. 1), the elastic relation still can be used to form an 
elastic compatibility equation: 


iAM, 1 
———? = - (AM, — AM.) —=pdpr (11) 
dr r 2 


Equipped with Eqs. (3) and (8) corresponding to the 
loading state, and Eqs. (3) and (11) corresponding to 
the unloading state, we are in a position to solve the 
clamped plate problems. Starting from the moment 
distribution of incipient plastic flow, we assume a cer- 
tain amount of AM, = AM, at the center of the plate 
with an increase in distributed load Ap. Eqs. (3) 
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and (8) are integrated by the digital computer along 
the plate radius. Each time a set of moment incre- 
ments is obtained for a certain radial position, it is 
tested to see if the plate element is in a loading state. 
The test is carried out by computing the value of AF 
(2M, — M,) AM, + (—M, + 2M,)AM,. If the value 
is positive, then the plate element is in a loading state. 
The integration process is resumed. If at a certain 
station (e.g., point B, Fig. 5), the computed value is 
negative, then the plate element at the station is in an 
unloading state. The transition point from loading 
to unloading lies somewhere between this station and 
the previous station. The location of the transition 
point is obtained by shortening the length of integra- 
tion starting from the previous station until the above- 
mentioned testing value equals zero. From this 
transition point on, the plate is in an unloading state 
Eqs. (3) and (11) are used accordingly. For the un- 
loading portion of the plate the largest value of the 
equivalent moment 


Mnaz = “M? — M,M, + M,? 
is recorded. In future tests loading can take place 
only when the value Mm: is exceeded. The integra- 
tion and test process continues until a station corre- 
sponding to a loading element is reached (point Q, 
Fig. 5). Then the transition point is found in a simi- 
lar manner, and Eqs. (3) and (8) are used again be- 
cause this portion of the plate is in a loading state. 
The integration process continues to the edge of the 
plate. 
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™ =1.50 
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Fic. 4. The circumferential bending moment distribution of the 
rigid, work-hardening plate. 











FIG 


ng 





PLASTIC BENDING OF A WORK-HARDENING CIRCULAR PLATE 819 


In the integration process, a set of computed A/,, 
AM, data is considered valid if condition (10) is satis- 
fied by the moment increments at the clamped edge. 
In order to reach a set of acceptable moment increment 
data with reasonable accuracy, the assumed A, value 
at the center of the plate is varied intentionally. In 
response to the variation in AM, value, some sets of 
(AM,, AM,) data at the clamped edge undershoot the 
limiting line (17, — 2M, = 0) which is plotted in Fig. 2, 
while some other sets of (AM,, AM,) data overshoot 
the limiting line. In case of undershooting, the dAM, + 
dr value as given by Eq. (8) is bounded so that there is 
no undue amount of error introduced. In case of over- 
shooting, the computed amount of overshooting is on 
the safe side because the extrapolation process is based 
essentially on the initial value of dAM,/dr along the 
integration step. The actual set of (AM/,, AM,) data 
lies somewhere between the two groups of extreme cases. 
When the gap between the two groups of data are 
closed, a final set of (AM,, AM,) data is reached which 
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Fic. 5. Possible moment curves of a work-hardening plate 


TABLE 1 
Radial and Circumferential Moments of a Clamped Circular 
Plate of Rigid, Work-Hardening Material Under Uniform 
Transverse Load 


(a) Radial Bending Moment M,/M, 

































































satisfies condition (10) to within a predictable amount pb/por/R = 0 0.2 0.4 0.6 0.8 1.0 
of error introduced by the numerical process. 1.00 1.000 0.935 0.725 0.336 —0.269 —1.155 
. 1.10 1.267 1.176 0.895 0.4380 —0.257 1.245 
1.20 1.488 1.3838 1.050 0.514 —0.254 —1.347 
1.30 1.684 1.569 1.192 0.590 —0.259 —1.455 
Numerical Data 1.40 1.862 1.737 1.322 0.657 —0.27 —1.571 
1.50 2.021 1.888 1.439 0.713 -—0.292 —1.698 
For a work-hardening plate, we integrate Eqs. (3) (b) Circumferential Bending Moment M,/M, er 
and (8) or Eqs. (3) and (11) by the digital computer a ; — 
, Y ; if p/por/R = 0 0.2 0.4 0.6 0.8 1.0 
to obtain the moment increments AM,, AM,, for 1.00 1.000 1.054 1.141 1.125 0.888 —0.577 
certain load increments Ap. After the moment incre- 1.10 1.267 1.262 1.261 1.213 0.889 —0.623 
Gita owas. ; a ee 1.20 1.488 1.462 1.385 1.299 0.935 —0.673 
ments and the enCTEMENTS ol moment derivativ es are 1.30 1.684 1.648 1.511 1.381 0.976 —0.727 
added to their initial values, the method is repeated for 1.40 1.862 1.820 1.636 1.461 1.010 —0.785 
: a : if 50 2.02 976 756 535 37 —0.8 
further increase in the uniform load Ap on the plate. = ia i i te oe 
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Fic. 6. 


The radial and circumferential curvature distribution of a clamped circular plate of rigid, work-hardening material under 


uniformly distributed transverse load. 
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Fic. 7. The deflection curves of the clamped circular plate 
of rigid, work-hardening material under uniformly distributed 
transverse load. 
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The results thus obtained are tabulated in Table 1 and 
plotted in Figs. 2—4. 

As can be seen from the given results, corresponding 
to successive load increments, nowhere does a plate 
element go through the unloading state. During the 
computation process, unloading of a plate element is 
encountered only when a wrong AM, value is assumed 
at the center of the plate. 

Substituting the obtained moment distribution for 
the work-hardening plate into flow rule (6), the plate 
curvature changes are calculated. The resulting data 
are given in Table 2 and are plotted in Fig. 6. The 
plate deflection is obtained by integrating r-k, along 
the radial axis. The deflection curves are plotted in 
Fig. 7. 


Conclusions 


The problem of a work-hardening plastic material is 
characterized by its nonlinear behavior. The plastic 
strain corresponding to a loading state is dependent on 
the stress history. In the unloading state, a quite 
different relation governs the strain change of the ma- 


terial. Thus, for a structure made of work-hardening 





TABLE 2 
Radial and Circumferential Curvatures of a Clamped Circular 
Plate of Rigid, Work-Hardening Material Under Uniform 
Transverse Load 





(a) Radial Curvature Ck,/M, 








p/por/R=0 0.2 0.4 0.6 0.8 1.0 
1.10 0.133 0.093 0.022 —0.014 -—0.029 —0.069 
1.20 0.244 0.184 0.053 -—0.024 -—0.057 —0.144 
1.30 0.342 0.269 0.091 -—0.031 —0.086 —0.225 
1.40 0.431 0.348 0.131 —0.036 —0.114 —0.312 
1.50 0.511 0.419 0.170 —0.040 -—0.140 —0.407 





b/por/R=0 0.2 0.4 0.6 08 1.0 
1.10 0.133 0.127 0.093 0.061 0.041 0.000 


(b) Circumferential Curvature Cky/M) 








1.20 0.244 0.238 0.184 0.124 0.083 0.000 
1.30 0.342 0.339 0.271 0.185 0.124 0.000 
1.40 0.4381 0.480 0.352 0.245 0.166 0.000 
1.50 0.511 0.512 0.428 0.302 0.205 0.000 
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material, the strains at various locations may respond 
according to different flow rules depending on the local 
stress conditions. It is this variation in response which 
makes the problem mathematically unwieldy. At 
the same time, a digital computer, with its ability to 
handle the computation according to local conditions, 
proves to be a convenient tool to handle this type of 
problem. 

For a problem where both loading and unloading take 
place in the work-hardening structure, interpolation of 
the stress (moment) values at the transition point is 
necessary. The errors induced by the interpolation 
can be limited to a predetermined magnitude by con- 
trolling the length of integration between two neigh- 
boring stations. 

To solve the differential equations by the digital 
computer, the Runge-Kutta method is used.’ The 
lengths of integration between two neighboring sta- 
tions are varied along the radius depending on local 
Eight significant figures are used for all 
the data. A truncation error in the order of 0.0005 MJ, 
is expected for both radial and circumferential moment 
This error is accumulated after more 


conditions. 


distribution. 
than 200 integrations are performed. 

While the compatibility equation (8) is partially 
based on the flow rule (6), the plate work-hardening 
constant C does not appear explicitly in the former equa- 
tion. This phenomenon indicates that, under the 
conditions specified in the paper, the moment dis- 
tribution is dependent on the load increments, the yield 
curve, and the physical dimensions of the work-hard- 
ening plate; it is independent of the degree of work- 
hardening of the plate material. Indeed, if we start 
from the state of incipient plastic flow and assume 
that the plate work-hardening constant C approaches 
zero, then the problem is reduced to that of a rigid plas- 
tic plate. As a result, in Figs. 6 and 7, the first sets 
of curves next to the curves corresponding to the state 
of incipient plastic flow (p/p) = 1), can be interpreted 
as the flow and deflection patterns of the rigid plastic 
plate with some modification of the coordinates. 

Fig. 4 indicates the circumferential moment dis- 
tribution of the work-hardening plate. Its nonlinear 
response to successive increasing loads is obvious. 
An examination of Figs. 3 and 4 shows that in some lo- 
cations the radial or circumferential moments do di- 
minish in magnitude with increasing load on the plate. 
The locations, where decreases in either radial or cir- 
cumferential moments are encountered, are staggered 
in such a manner that the plate elements stay in the 
loading state. 
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Inviscid Flow With Nonequilibrium Molecular 
Dissociation for Pressure Distributions 
Encountered in Hypersonic Flight’ 


MARTIN H. BLOOM* ann MARTIN H. STEIGER** 
Polytechnic Institute of Brooklyn 


Summary 


One-dimensional inviscid nonequilibrium flows of a_ two- 
component model gas are studied for prescribed pressure vari- 
ations and an average reaction rate based on recent data for 
oxygen recombination. These flows are interpreted in relation 
to the flow along streamlines around blunt hypersonic bodies. 
Assuming equilibrium conditions in the subsonic region, it is 
estimated that the flow in the initial supersonic expansion region, 
which is approximately of Prandtl-Meyer character, will be 
chemically frozen with respect to the molecular dissociation of 
the primary components under the hypersonic, high-altitude 
flight conditions considered. The flight conditions consist of 
flight velocities between 15,000 and 25,000 ft/sec at altitudes 
between 154,000 and 246,000 ft. Furthermore, on bodies of small 
surface inclination beyond the nose, the flow will continue to be 
effectively frozen for at least 20 ft downstream of the nose. These 
conclusions may lead to the simplification of procedures for theo- 
retical calculation and testing. 

The problem of distinguishing a dimensionless length-reaction 
rate parameter, which characterizes the extent of departures 
from equilibrium or from frozen behavior in the flow fields of 


interest here, is discussed. 


Symbols 


(Units shown are only typical; consistent equivalent units may 
be used in calculation. ) 
mass fraction of species 7 


Gc = 

C = k{1 + a@)pa/m,?, parameter defined in Eq. (14) 

d = dissociation energy per molecule (ergs/molecule ) 

Dy) = dissociation energy per mole (ergs/mole ) 

D = dissociation energy per unit mass of molecules (ergs 
gm or cm?/sec?) 

D = D/(u.2/2) 

Dy = diffusion coefficient of binary mixture 

e = internal energy per unit mass 

f = partition function 

1 = enthalpy per unit mass 

i = ¢%,,*/2) 

% = #{#,*/2) 

k = Boltzmann’s constant, 8.62 X 10-> ev/molecule °K 
or 1.380 X 10~' ergs/molecule °K 

k, = recombination rate parameter (cc/mole )?/sec 

L = characteristic length of boundary 

m, = atomic mass (gm/mole) 

m, = molecular mass (gm/mole) 
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molar concentration per unit volume (moles/cc) 


n = 
p = pressure 
D = bf PotUe?/2) 
r = frontal radius of flat-nosed body 
Yn = nose radius or typical dimension 
R = Ro/2m, gas constant for undissociated diatomic gas 
Ry = 1.987 cal/mole °K, universal molar gas constant 
5 = distance measured along a streamline 
s = §/L,nondimensional distance 
t = time 
T = temperature, T = RT/(u,2/2) 
u = velocity along streamline 
u. = free-stream (flight) velocity 
V = local mean velocity vector 
V; = diffusion velocity of species i (Eq. (7)] 
w; = rate of production of species 7 (moles/cc sec) 
a = mass fraction of atoms (identical with c,) 
y = ratio of specific heats 
p = mass density 
p = p/P 
pa = characteristic density for a molecular species [Eq. (6)] 
y = streamline deflection angle 
¢g¢ = (LCp../taopa)™ 
Ga = (LCpq?/tapa)™ 
Subscripts 
i = property of a member of species 7 
1 = atomic species, or conditions at s = 1 
2 = molecular species 
s = stagnation values 
0 = see complete symbol 
00 = undisturbed flight conditions 
e = equilibrium conditions 


aor sonic = sonic conditions 


horiz = conditions on surface parallel to uv, 


Preliminary Discussion 


| ¥ IS USEFUL TO sTuDY the effect of nonequilibrium 
molecular dissociation on flow variables, such as the 
dissociation fraction, density and temperature, when 
the flows are subjected to prescribed pressure distribu- 
tions and when transport properties are negligible. 
In particular, these parameters are of interest when 
the pressure distributions correspond to those existing 
over blunt bodies in hypersonic flight. 

For reacting gases of this type, complete flow field 
calculations giving the pressure as well as the other 
properties over blunt hypersonic bodies require com- 
plex procedures in both the subsonic and supersonic 
regimes. The present procedure is much more amena- 
ble to a parametric study covering a wide range of 
flight conditions. 
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Pertinent Literature 


Recent representative basic literature concerning 
chemically-reacting gas flows includes the work of 
Penner,' Hirschfelder, Curtiss and Bird,’ the compila- 
tion edited by Lewis, Pease and Taylor,’ and a brief 
Although these are oriented 
combustion 


review by Emmons.‘ 
somewhat toward the consideration of 
processes, the general concepts—consisting of a syn- 
thesis of the conservation laws of flow, thermodynam- 
ics, and chemical kinetic principles—apply to the 
present problem as well. 

More directly relevant to the present problem in- 
volving dissociative reactions are the works of Chu,5 
Resler,® Lick,” * Clarke,? and Li,’° which review the 
governing equations, examine aspects of wave propa- 
gation (particularly the role of the “‘frozen’”’ and equilib- 
rium sound speeds), and discuss the solution of super- 
sonic flow fields by the method of characteristics. 
Lick” * moreover has presented an adaptation of the 
inverse method for computing the subsonic flow behind 
a prescribed detached shock and has calculated some 
typical hypersonic flow properties. Wood and Parker"! 
have studied the structure of a centered rarefaction 
wave subject to a single internal relaxation process 
(such as molecular vibration) but without chemical 
reaction (such as dissociation), utilizing the method of 
characteristics. On the other hand, the influences of 
dissociation reactions on simple shear and boundary 
layer flows, including transport effects, have been 
treated by Broadwell,'? Clarke'*® and Jarre.'4 Napoli- 
tano’’ has examined a centered rarefaction flow with 
nonequilibrium dissociation by means of expansion pro- 
cedures in terms of the radial coordinate, assuming 
chemically-frozen flow near the origin and equilibrium 
flow far from the origin. 

Since the approach of this paper involves prescribing 
a pressure distribution along a streamline or streamtube 
and computing the consequent remaining flow variables, 
it is essentially a one-dimensional analysis. Calcula- 
tions of one-dimensional nature have been made by 
Wood,'* Evans,” Freeman'* and Bray'® for the flow 
behind normal shocks, for which the streamtube area 
remains constant. One-dimensional calculations have 
also been made by Heims,” Bray,'* Hall?! and Browne”? 
for the flow in nozzles with prescribed area distributions. 
Heims” moreover interpreted the nozzle as corre- 
sponding to a streamtube along a blunt body, but sub- 
sequently assumed a rather arbitrary converging- 
diverging streamtube geometry for a single numerical 
example, and did not carry the analogy further; a ma- 
jor part of his paper was devoted to a survey of 
reaction rates for oxygen dissociation. Rudin® has 
also studied this problem. A comprehensive review 
paper concerning chemical reactions in nozzle flow has 
been composed by Penner and Jacobs.”? 

It may be remarked that Wegener,** in an experi- 
mental study of the dissociation and recombination of 
N2O, in a working fluid of nitrogen in a nozzle, has re- 
ported agreement with Bray’s method of predicting 
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deviations from thermodynamic equilibrium and com- 
ponent freezing; the thermodynamic and reaction rate 
data for N:O,; were known quite accurately for this 
study. 

Finally, we may cite Freeman’s'* analysis of the non- 
equilibrium flow field in the low Mach number region 
of a blunt body of revolution. An essential simplifi- 
cation in his analysis resulted from the determination 
of the pressure field explicitly, prior to the execution 
of the remainder of the calculation. Thus the proce- 
dure was reduced ultimately to the calculation of prop- 
erties along a streamline, and a single set of calculations 
was carried out for a given free-stream density and 
various values of the reaction rate parameter. The 
pressure distributions within the approximation used 
by Freeman were given by the ‘‘Newtonian-plus- 
centrifugal force’ value, which was valid over a sphere 
only up to points a little beyond the sonic region (that 
is, for } < 60°, where # is the angle between the surface 
normal and the body axis). At # ~ 60° the pressure 
diminished to zero due to the centrifugal effect accord- 
ing to the theory used. This unsatisfactory predic- 
tion for the pressure has recently been eliminated by 
Freeman” through an improvement in the theory; 
the revised pressures are in agreement with the ‘‘modi- 
fied Newtonian’”’ which was previously 


strictly empirical. 


prediction 


A Type of Similitude 


It is important to note that in the aforementioned 
inviscid analyses, the reaction rate constant, that is, 
the numerical factor which is directly proportional to 
the reaction rate and may represent a characteristic 
time, appears only where it may form a product with 
the geometric coordinates or with a reference length 
of the boundary. Therefore only the product of the 
rate constant and the reference length need be specified 
in a given calculation. The rate constant may be 
arranged to appear in the differential equations or in 
the boundary conditions, as desired. 

This similitude is somewhat restrictive for practical 
application, since if one wishes to obtain the flow fields 
about a series of geometrically similar bodies for a given 
value of the rate constant, a separate calculation must 
be made for each scale length. Likewise, if a single 
body is considered for a set of different rate constants, 
a separate calculation must be made for each rate con- 
stant. On the other hand, use of the above similitude 
permits us to determine, from the flow about one body 
with a given rate constant, the flow with a second rate 
constant around one other geometrically similar body. 

In the blunt-body calculations of Freeman" the 
rate constant is retained in the differential equations, 
whereas in the calculations of flow behind normal 
shocks, the flow is considered to be one-dimensional 
with constant cross-sectional area and a geometric 
boundary condition is not required; hence the rate 
constant can be removed from the differential equa- 
tions and finally appears as a scale factor of the stream- 
wise coordinate. 
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INVISCID 


It is suggested that more general similitude relations 
may be obtainable for certain special classes of flow, 
such as those amenable to small-perturbation treat- 
ment (e.g., linearized supersonic flow or some hyper- 
This possibility will not be considered 
Studies such as those of Moore 


sonic flows). 
further in this paper. 
and Gibson* and Vincenti,*® who treated the linearized 
nonequilibrium flow over a wavy wall, are of interest 


in this regard. 


Choice of Pressure Distributions 


For the purposes of this paper, it is considered rea- 
sonable to assume that the pressure distributions over 
blunt bodies of simple type under hypersonic flight 
conditions will not be altered severely by nonequilibrium 
thermodynamic processes or indeed by real gas effects 
in general, if an effective constant value of the specific 
heat ratio or adiabatic exponent y is used in calcu- 
lating the pressure field, at least as a first approxima- 
tion. Thus pressure distributions selected on the 
basis of thermodynamic equilibrium, frozen composi- 
tion, or estimates neglecting real gas effects, for ex- 
ample, can be used as a basis for evaluating all other 
properties over the body. Viewed in another way, the 
prescribed pressures can be considered to correspond 
in actuality to body shapes which are somewhat, but 
not greatly, different from those which originally fur- 
nish the pressures. 

Correspondingly, 
streamlines can be prescribed according to conditions 
behind the bow shocks of hypersonic bodies. In this 
paper only stagnation conditions close to those existing 
behind normal shocks are considered in the calculations. 


stagnation conditions along the 


Thermodynamic Functions and Reaction Rates 


An idealized gas is selected here to represent the ap- 
proximate behavior of real air. The gas is considered 
to be composed of only a single type of symmetric 
diatomic molecule which is free to dissociate into two 
atoms as a result of two-body collisions. The atoms, 
on the other hand, may recombine into molecules as a 
result of three-body collisions. The vibrational energy 
of the molecules is taken into accounr in the thermo- 
dynamic relations in an average fashion after Lighthill* 
who suggested the inclusion of one half the classical 
vibrational energy for all molecules, with a consequent 
simplification of the thermodynamic relations. This 
simplified thermodynamic analysis as well as more com- 
plete ones have been discussed in detail in most of the 
previously cited references (e.g., 5, 7, 9, 15, 27) and 
need not be repeated here. However, it can be stated 
that Lighthill’s analysis was made under conditions of 
thermodynamic equilibrium. For nonequilibrium dis- 
sociation it is still assumed that each subsystem, one 
consisting of undissociated molecules and the other of 
atoms, is in equilibrium in itself, and the thermody- 
namic functions are calculated on this basis. 

The approximation of real air by means of a simple 
diatomic gas is not too farfetched, since nitrogen disso- 
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ciation will be dominant in the range of state conditions 
to be considered here, while the molecular weights of 
oxygen (O. = 32) and nitrogen (N2 = 28) are quite 
close to each other. Furthermore, in the calculations 
to be made here it will be possible to make an empirical 
correction for combined nitrogen and oxygen dissocia- 
tion effects by estimating an ‘effective dissociation 
energy’ for each calculation instead of assuming the 
reaction to be governed by the dissociation energy 
of nitrogen or oxygen alone. This effective value will 
be determined by requiring that the simplified mass 
action law of the model gas be satisfied for the true 
equilibrium state values of air given, say, by one of the 
recent tabulations or by Mollier charts for a reference 
equilibrium state. This effective 
values somewhat higher than the dissociation energy 
of oxygen, due to the contribution of the higher dissoci- 
ation energy of nitrogen (Table 1). The procedure 
will be discussed further in the section dealing with 
the calculations. It appears also that the reaction 
rate values for oxygen and nitrogen are not substantially 
different from each other (see Zinman*®). Therefore 
the oxygen values used here are assumed to be reason- 


procedure yields 


ably valid for air reactions. 

The influence of ionization and radiation phenomena 
on the thermodynamic behavior of the primary air 
components will be neglected. These effects are not 
expected to be large. 

The pertinent thermodynamic relations are summa- 
rized as follows: 


Equation of state: mixture of symmetric diatomic 


molecules (m. gm/mole) and their atoms (m, = m»/2, 
gm/mole) treated as perfect gases. 
p= p(1 + a)RT (1) 
where 
R = Ri/m, gas constant per unit mass of molecules 
(cal/gm °K) 
Ry = kNo, gas constant per mole (1.987 cal/°K 
mole) 
k = Boltzmann constant, gas constant per mole- 


cule (cal /molecule °K) 
N, = Avogadro constant (molecules/mole) 
mass fraction of atoms in mixture 


! 
| 


a 
Internal energy per unit of mixture, e, including one 
half the vibrational energy of the molecular species 
(the Lighthill idealization appears in this relation): 


e = 3RT + Da (2) 
where 
D = D,/m:, dissociation energy per unit mass of 
molecules (cal/gm) 
Dy = dNp, dissociation energy per mole of molecules 
(cal/mole) 
d = dissociation energy per molecule (cal/mole- 


cule) 


Typical values of the dissociation energy are listed in 
Table 1. 
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TABLE 1 
Typical Values of Dissociation Energy 





Quantity Units Nitrogen (Ne) Oxygen (O2) 
a ev/molecule 9.76 5.08 
ergs/molecule 15.65 X 10712 8.14 xX 107" 
Do = Ned _ ergs/gm mole 9.45 x 10% 4.91 x 10% 
D =D)/m_ ergs/gm 2.95 < 10" 1.536 <* 101! 
dyne-cm/gm 
(gm-cm?/sec?) / 
gm 
(cm/sec )? 
ft-lb/slug 3.18 xX 108 1.651 * 108 
ft?/sec? 
d/k _.- wees 113,000  _—_—59, 000 


Enthalpy per unit mass of mixture, 7: 
i=e+ (p/p) = (4+ a)RT + Da (3) 


From the standpoint of a chemical reaction whose 
reaction rate must be determined, the dissociation and 
recombination of an idealized symmetric diatomic 
molecule may be treated as a simple binary reaction of 
the form (e.g., see references 9, 12, 18) 


k, 
Ag + As=— 2414+ Az (4) 
k, 


where A» represents a molecule, 4; an atom, and A; 
either an atom or molecule which participates in a 
collision and emerges unchanged in form. The chemi- 
cal kinetics of such reactions is well known (e.g., refer- 
ences 9, 12, 13, 18, 19) and readily yields the functional 
form of the rate relation as follows: 


a E + a) | ? | \ ,—D/RT _ ] 
W, = | - : Pa (1 —- ape ——a 
p m," Pa Pa 
(5) 
where 
w, = net rate of production of atoms, 
(moles/ce sec) 
a = mass fraction of atoms in the mixture 
m, = gram atomic mass (gm/mole) = m2/2 
mM, = gram molecular mass (gm/mole) = 2m 
pa = (m;/2)(fi?/fe)(gm/cc) 
fife = partition functions of atoms and molecules, 
respectively 
k, = recombination rate function (cc/mole)?/sec 


Eq. (5) is still general in the respect that no approxi- 
mation has yet been made concerning the nature of the 
partition functions f; and fs contained in pg. The 
Lighthill”’ idealization, taking into account one half 
the classical vibrational energy, is introduced if vari- 
ations in pa are neglected, with the result that fi? ~ fo. 
The assumption pg = constant will be used henceforth 
in this paper. The following average values are cited 
from the work of Lighthill: 


. = $e J 
for oxygen: pa = 150 gm/ccl jog & F< 7000°K 
for nitrogen: pa = 130 gm/cc 
(6) 
. t / 
Sercuygen: fa = 180 gmi/ect 5000 < T < 7000°K 
for nitrogen: pa = 125 gm/cc! 
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The latter values will be used in the calculations to be 
made here. 

It will be necessary to combine Eq. (5) with the 
equation of mass conservation of atomic species for use 
in the present study. The mass conservation equation 
for each species can be written as: 


On, > > oy 
+ V-[n(V + V,)] = w, (7) 
Ot 
where 
V = local mean velocity 
Vi = — (Dw/ci)Vc;, random or diffusion velocity of 


species 7 (to be neglected here) 
Dy = diffusion coefficient of a binary mixture 


c; = :M;i/p, mass fraction of species 7 (note: 
qQ= a) 
w; = rate of production of species 7 (moles/cc sec) 


By taking into account the overall continuity equa- 
tion formed from the set of Eq. (7), namely, 
l Dp i> 


Dr +V-V=0 (8) 
p 


Eq. (7) can be arranged as follows: 
De, mW, 


- DyeV c:) = (9) 
Di , (pD, 


Neglecting the diffusion term and considering the 
atomic species we obtain from Eqs. (5) and (9) 


Da = E = a) ee £ t — ae . E a | 
Dt m2 Pa Pa 


(10) 


The Recombination Rate Parameter k, 


An extremely important factor in this investigation 
is the actual magnitude of the parameter k, which 
finally establishes the reaction rate. Experimental 
data concerning oxygen dissociation recently has been 
presented by Byron** and Matthews.*’ These data are 
believed to be comparatively reliable; the two sets of 
results are in essential agreement with each other and 
probably establish the value and temperature depend- 
ence of k, within multiplying factors of about 3. The 


99 


following relation is given by Matthews:* 
k, = 8.4 X 10'(7/3,500°K)—*(cc/mole)?/sec (11) 


while Hall*! cites the following value based on the work 
of Byron and of Matthews: 


k, = 1.2 X 10% (7/3,000°K)—* X 


(cc/mole)?/sec (1.0 < s< 2.5) (12) 


Reviews of some values of k, obtained or estimated 
prior to the above recent work have been given, for 
example, by Heims,”® Bray! and Clarke.'* The vari- 
ous values cited differed from each other by multi- 
plying factors as high as 1,000 and will not be considered 
further here. 

According to the previously mentioned recent study 
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by Zinman,” the recombination rate parameter for 
nitrogen is close in value of that of oxygen. 

For the present study we shall use the value of k, given 
in Eq. (11), but will neglect the temperature dependence 
effect. Thus with 


k, = 8.4 X 10!4 (cc/mole)?/sec 

Pa = 130 gm/ce . 
: (13) 

m = 16 gm/mole 

lt+anrl4 


we obtain the following estimate of the leading factor 
in Eq. (10): 


C= k,(1 + a)pa/m,? = 6 X 10'4(cc/gm)/sec (14) 
= 3 X 10'4(ft?/slug)/sec 


t must be remembered that the cited values of k, 
were obtained in pure oxygen or oxygen-argon mixtures. 
They may differ somewhat for oxygen dissociated in 
air, but are the best values available at present and may 
be expected to yield satisfactory engineering estimates. 
The value of C given in Eq. (14) is used for the present 
calculations. 


System of Equations 


The following equations govern the nonequilibrium 
steady flow along streamlines in an idealized gas whose 
properties resemble those of oxygen as discussed in the 
previous section: 


du ldp . 
momentum: u = 0 (15) 

ds pds 

me 
energy: t+ — = (16) 
continuity: puA = const. (17) 
state: pb = pRT(1 + a) (18) 
enthalpy: t= (4+ a)RT + Da (19) 


reaction rate plus atom conservation: 


d " 
u = = LCp c — aje~P/RT . | (20) 
ds Pa 


where s=3/L 

Note that the continuity relation may be used in a 
final step to yield the cross-sectional area distribution 
A(s) of the streamtube. The momentum and energy 
equations retain their usual inviscid, adiabatic form. 
Also, the product LC appears in Eq. (20) as previously 
discussed. 

Since hypersonic flows only are being considered here 
we may write 


i? = u.*/2 (21) 


By introducing the dimensionless parameters defined 
in the notation list, we may rearrange Eqs. (15)—(20) to 
yield the following working forms: 
4+ a)p - 
= et. + Da (22) 
(l+ a)6 
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Fic. 1. Fractional change in a at sonic and fully expanded 
regions 


T = RT/(u 7/2) = p/p(l + a) (23) 
d(a*) | (1+ a)(1 — a — Da) 1 dp 
‘4 = (24) 
ds (4 + a) pb ds 
_ da (A *) 
u aod, lee p X 
ds i 
Sy ) | Mat © | Ps a} ) 
- exp] — _ pa? 25 
\ ” I 1 — a? — Da Pa = ~ 
where it is noted that D = D/(u.?/2) and thus de- 


pends on flight conditions. 


Extent of Variation in a for Equilibrium Flow 


To assess the range of variations in the mass fraction 
of atomic species, a, that could be expected over a 
range of altitudes and velocities, a series of calculations 
was made assuming chemical equilibrium, and conse- 
quently isentropic expansion, from various stagnation 
conditions behind a normal shock corresponding to 
flight velocities from 15,000 to 27.000 ft/sec and alti- 
tudes from 140,000 to 300,000 ft. Shown in Fig. 1 
is the percentage change from the stagnation value 
a, for an expansion to the sonic pressure and for a fur- 
ther expansion to P/Psonie = 0.1142, a reasonable pres- 
sure ratio for expansion to horizontal (surface de- 
flection zero with respect to u..) on a blunt body. The 
values shown have been obtained for real air. 

It is interesting to note that the changes in a are 
limited rather narrowly to the following ranges over 
the entire altitude-velocity regime considered: 


(Qs — Qonie)/A ~ 0.05 to rp | 
(Qs — Qnoriz)/A ~ 0.22 to 0.29 (26) 
(Qonte — %nortz)/A@ ~ 0.17 to 0.24 
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Fic. 2. 


These changes can be interpreted as representing 
upper bounds for those which will occur in nonequilib- 
rium flow. At the other extreme of flow under condi- 
tions of chemically-frozen composition, the correspond- 
ing change is zero. 


Selected Pressure Distributions 


‘‘Prandtl-Meyer’’ Variation 

To achieve a desirable degree of generality in the 
parametric study, the flow region downstream of the 
sonic region on blunt bodies will be examined (Fig. 12). 
In this part of the flow sharp decreases in pressure 
occur together with velocities which are substantially 
larger than those in the initial subsonic region. There- 
fore relatively marked deviations from thermodynamic 
equilibrium may be expected to appear. 

It is considered reasonable to expect that the pres- 
sure along a streamline in this regime will resemble 
that given by a conventional Prandtl-Meyer expan- 
sion for constant y. Of course this is necessarily a 
rather rough approximation, neglecting the detailed 
features of the complicated flow field in this region. 
However it is considered suitable for the demonstration 


sought here. It will hold most closely for streamlines 


over cylindrical or spherical noses, particularly those 
which emerge from the shock in the subsonic region 
and subsequently remain close to the body. 


The 
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length over which the assumed pressure variation is 
assumed to be distributed on these streamlines is of the 
order of one nose radius starting from the sonic point 
or each streamline; that is, for spherical roses, 


L= O(rn) (27) 


Of course the streamline lengths within the shock layer 
some distance away from the body will be somewhat 
longer than those near the body, due to the inclination 
of the sonic line. These increased lengths, which will 
not usually exceed the value 1.5 7,, approximately, still 
fall within the order designated by Eq. (27) for present 
purposes. 

In the case of flat-nosed bodies with sharp rims or 
small rim radii, the streamline lengths in the initial 
supersonic expansion region of interest here are ex- 
pected to be on the order of the shock-detachment 
distance or less, being included within an ‘expansion 
fan’’ angle of about one half radian. The shock- 
detachment distance, which actually may be correlated 
on the basis of the shock radius, has been shown™ for a 
flat-nosed body of revolution to be roughly equal to 
one half the body radius. Therefore for this type of 
body at least, the maximum characteristic flow lengths 
are roughly equal to 1/4 the body radius, or, for flat- 


nosed bodies: 
L = 0/4) (28) 


Of course, streamlines directly adjacent to the sharp 
rim will pass through the expansion fan in an infini- 


tesimallength. In this case we specify 


L «QO (29) 


The pressure variation through a Prandtl-Meyer 
expansion is plotted in Fig. 2 versus streamline deflec- 
tion angle v for y = 1.4 and y = 1.1, the latter value 
being representative of hypersonic conditions. For a 


surface with constant radius of curvature r,, in this re- 


gion, the angle v equals the arc length s = $/r,. The 
following expression fits the indicated variation for 
y = 1.1, and is used for the calculations made here: 

p 


= ().2284 + 1.7716(1 — s)?-4 (30) 
(Deonic/ Ps) 


Constant Pressure 


A type of pressure distribution which can be con- 
sidered in a relatively simple fashion is that of 


p = const. (31) 


for which u and / are also constant, and only one differ- 
ential equation [Eq. (25)] need be integrated for a. 
Constant pressure conditions may be encountered 
over certain types of afterbody such as the conical skirts 
of blunted cones. In the case of a zero cone angle 
(blunted cylinder) the surface pressure generally de- 
creases, especially near the nose junction; hence the 
body must be flared somewhat to produce constant 


pressure. 
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“Exponential’’ Variation 


Another special pressure distribution fulfills the 
condition : 
d(In p) /d(a?) A = const. (32) 
or 
p = Be**, B = const. (33) 
for which Eq. (24) becomes 
l+a a 
+ ( Jo — a7 — Da)A = 0 (34) 
f+a 


where 4 and B are evaluated at reference conditions of 


p,aiand a. Equation (34) shows that 4 < 0. 

In this case again only Eq. (25) need be integrated 
for a. Equation (33) has been shown in typical calcu- 
lations* to yield decreasing pressures along afterbodies 
such as those of blunted cones of small cone angle or of 
blunted cylinders. 

Pressure distributions of the two latter types can 
be used to examine the nature of the return to equilib- 
rium of flows which have deviated from equilibrium 
in the nose region. The calculations are convenient 
particularly if automatic computing equipment is un- 


available. 


Calculations and Results 


“Prandtl-Meyer’’ Pressure Distribution 


Integration of (24) and (25) with the pressure dis- 
tribution of (30) yields the flow properties along a 


* Not presented in this paper 
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streamline. In the present calculations the following 


values are used in (25): 


pa = 250 slugs/ft® { 
C = 3 X 10" ft*®/slug secf 


The initial conditions for the calculations are esti 
mated to resemble conditions in the sonic region of a 
blunt body, assuming the flow prior to the sonic re- 
gion to be in equilibrium. The sonic conditions are 
related to flight conditions by utilizing a Mollier chart 
for real air in equilibrium, assuming an isentropic ex- 
pansion from stagnation conditions behind a normal 
shock. The values utilized are listed in Table 2. 

For each initial (sonic) state it is desired to have the 
equilibrium mass action law of the model gas exactly 
This 


is necessary for the integration procedure to be de- 


satisfied at the prescribed values of p, @ and 7. 
scribed. The required expression is obtained by set 
ting the term in braces in (25) equal to zero and ex- 
pressing it in the form: 

D(A + a) 
i — Da 


(| *) pa- 

= () 
l + a ; = Da 
For prescribed initial values of ), a and 7, a value of 
This value is termed the 


pall — a) exp] — 


(36) 


D is obtained from (36). 
“effective value of D,’’ denoted by D. +5. 
sures that (36) is satisfied at the initial point of the 
calculations. Table 
2 along with the corresponding dimensional values of 
D.5;, which are seen to be generally higher than the 
dissociation energy of oxygen and closer to that of 


Its use en- 


Values thus obtained are listed in 


nitrogen (Table 1). 


TABLE 2 
Flow Conditions for ‘‘Prandtl-Meyer’’ Expansion 


Stagnation Conditions 











Altitude— ’ ——Normal Shock 

Po» (42° (z : he we 

Alt., slug, he pa), Pal Po ft? slugs 

Reference Gus ft, Hu. Pe 1/ft, U « Pa Ps sec?, Ts, ft3, 
number ft/sec < 10-3 x 108 psfa «K 105 x 1073 xX 108 psfa x 1078 kK a x 105 
1 15,000 154 287 2.52 65.9000 57.40 1.1480 645.75 1.1556 4,962 0.232 3.3398 
2 200 54.0 0.423 2.3328 0.80 0.2160 121.50 1.1504 4,717 0.242 0.6645 
3 246 8.40 0.0562 0.0564 1.68 0.0336 18.900 1.1463 4,469 0.251 0.1099 
} 20 , 000 154 287 2.52 49 4214 43.05 1.1480 1,148.00 2.0306 6,235 0.419 4.1367 
§ 200 54.0 0.423 1.7496 8.10 0.2160 216.00 2.0254 5,782 0.438 0.8393 
6 246 8.40 0.0562 0.0423 1.26 0.0336 33.60 2.0213 5,340 0.456 0.1399 
7 25,000 154 287 2.52 39.5371 34.44 1.1480 1,793.75 3. 1556 7,159 0.679 4.8361 
8 200 54.0 0.423 1.4000 6.48 0.2160 337 . 50 3.1504 6,569 0.706 0.9842 
9 246 8.40 0.0562 0.0339 1.01 0.0336 52.500 3.1463 6,013 0.732 0.1624 

— — Sonic Conditions 

Alt - lay ? Pa, D ’ 
Reference ea ft, Be, ft?/sec?, m slug/ft' ies ™ ft?/sec? 
number ft/sec x 10 psfa x 10-* °K ; Qa XK ft/sec Dy, x 1078 
es ~ 15,000 154 351.4 1.052 4,584 0.219 2.0309 4,560 2.0559 2.3123 
2 200 69.5 1.052 4,411 0.229 0.4057 4,447 2.0922 2.3537 
3 246 10.7 1.053 4,217 0.238 0.0674 4,325 2.1362 2.4032 
4 20,000 154 666.4 1.884 5,959 0.391 2.5629 5,414 1.4391 2.8782 
5 200 124.1 1.888 5,543 0.410 0.4996 5,234 1.4465 2.8930 
6 246 18.9 1.894 5,150 0.429 0.0824 5,051 1.4523 2.9046 
7 25,000 154 1023.5 2.954 6,822 0.642 2.9376 6,347 0.99868 3.1209 
8 200 193.1 2.963 6,290 0.669 0.5920 6,118 0.99939 3.1231 
9 246 20.4 > 973 5,785 0.696 0.0999 5,887 0.99861 3.1207 
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Fic. 3. Mass fraction of atoms for equilibrium expansion of 
model gas and real air 


The values of p and 7 calculated from the model gas 
relations (22) and (23) for the aforementioned pre- 
scribed values of p, a, 7 and D4 do not coincide pre- 
cisely with the value for air in Table 2. This cannot 
be avoided in a comparison between air and the model 
gas, although either p or 7’ may be made equal to the 
air value if an effective value of the gas constant RK, 
denoted by R,,,;, is defined to satisfy either (23) or 
(24). This refinement is not utilized here. The 
initial values of p and 7 obtained for air are compared 
with those based on the model gas in Table 3. 

If either p or 7 had been matched to the air-equilib- 
rium values in satisfying (36), values of D,,, differing 
on the order of 10% from those utilized here would have 
been obtained. These would not be expected to alter 
the trends indicated by the present results. 

A Runge-Kutta integration technique on a Bendix 
G-15 digitai computer has been used for the calcula- 
tions. The initial conditions are such as to make da + 
ds = Oats = Oin (25). However, this zero slope is not 
discernible in the practical scale of the intervals em- 
ployed here to present the results. 

The question arises as to what the initial slope would 
be in the limiting case tending toward equilibrium as 
L— © or C— o in (25). One would expect the slope 
to approach that of the equilibrium flow, that is, to 
approach the value of da,/ds at s = 0; the latter value 
is not zero. The treatment of this limiting case repre- 


~ TABLE 3 


Comparison of Density and Temperature at Sonic Conditions for 
Air and the Model Gas 


“p (slugs /ft®) 








Air Model gas — r CR) — 
Condition * x 105 x 105 Air Model gas 
1 2.0209 2.3640 8,251 7,074 
2 0.4057 0.4910 7,940 6,681 
3 0.0674 0.0797 7,591 6,291 
+4 2.5629 2.8900 10,726 9,616 
5 0.4996 0.5730 9,977 8,909 
6 0.0824 0.0949 9,270 8,201 
a 2.9376 3.1460 12,280 11,492 
8 0.5920 0.6370 11,322 10,531 
9 0.0999 0.1080 10,413 9,624 





* See Table 2 for corresponding velocity-altitude values. 
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sents a singular perturbation problem in which ¢ = 
(LCp../u pa)! — O and da/sa remains finite; it is 
treated in some detail in a separate report.*' It is 
sufficient to state here that the results of reference 31 
indicate that for practical purposes it is still sufficient 
to formally permit the initial slope (da/ds),.9 to be 
zero even for the approach to equilibrinm in the present 
calculations. 

We must be careful to note here that the criterion for 
judging what is meant by ¢ — 0 for equilibrium (or, for 
that matter, what is meant by 1/¢ — 0 for frozen flow) 
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Fic. 4. Distributions of @ for ‘‘Prandtl-Meyer’’ expansion for 
equilibrium and nonequilibrium flow. 











for 
for 
yw) 


537 


37 


or 
















































INVIS 

70 l — 

! 
it — 

, : . — + 
ail 200000 Ft '° 0.00014 

| 25000 FT " 

100' O0014 

| | ~ 
) a re i ae 

| P y. 

P N 
088+—— “~ ——+ 
/ | 
A 
AIR EQUILIBRIUM 
| | | MODEL GAS 
0.50 1 +. EQUILIBRIUM 
(¢) 

— a 200,000FT! —Ssi 

—o . | | 

| il — 

' coer names (FN $ 
040 1100' OOOI75 
Q35;— 
a30r— /+—_+___{__ 

| AIR EQUILIBRIUM | 
| | MODEL GAS 
| | , EQUILIBRIUM 
0.25- | { — 
(e) 
TT ee 
15,000 Ft/ 
SEC 
1 
o25+- t — + — 4H — 
| | . é ) 
z 7 100' 0002333 


1000° 0.02333 





020+— — 
PP sina 
—~——_ a 
| | | 
“a 
|< MODEL GAS 
| | <AIR EQUILIBRIUM | ama senies 
| | | 
ow L—_—_ 1 1 ‘in 
fe) 0.2 04 0.6 08 LO 
Ss (d) 
Fic. 4. Continued. 


GT 


CID 


MT 

















































































































FLOW 829 
Q75 
246,000 FT | 
25,000 FY 
\ SEC 
0.70r | | 
a | 
RS 
0.65 WN 
A 
AN 
a6o ‘ 
‘a 
i 
Qs5+ . al 
AIR EQUILIBRIUM MODEL GAS 
EQUILIBRIUM 
0.50 
(i) 
0.50 
246,000 FT 
20,000 FT 
SEC 
0.45 
A 
0.40 = 
NS 
S N 
0.35 | = 
4 | - 
“AIR EQUILIBRIUM a 
| \ | MODEL GAS 
030 EQUILIBRIUM 
(h) 
0.30 T “ 
246,000 FT 
15,000 FT/ 
AEC 
0.25 
| — 
} — 
— 
tea = 
“ it - 
Q15 +— “air EQuILiBRiuM 
MODEL GAS 
| EQUILIBRIUM 
! 
0.104 1 - 
0 0.2 0.4 0.6 08 Oo 











(g) 
Ss 
Fic. 4. Concluded. 





830 JOURNAL OF THE AEROSPACE 










wee te oe 


L*i000'and EQUIL 
L=Lo* 





| | 
SL 


U.7 15,000 FT 
; | /SEC 








— 





“L# 100) 





| ___- 200,000 Ft | = 




















I. 
fe) T Us *20,000 Ft 
| | | EC. 
| 
od T — ————— ee 2 
2.0 | a | L=10’ 
Leo’ | 
| 
LO act 154,000 FT __ | 
u, *25,000 a 
SEC. 
| | 
ae ae | 
fe) a2 0.4 a6 os 1.0 12 


Ss 


Fic. 5 Velocity distribution for equilibrium and nonequi- 
librium flow (model gas) 


has not been established as yet. It will be seen that 
values of ¢ as high as 10 may yield flows which are 
virtually in equilibrium, whereas values of 1/¢ as low 
as 0.005 result in flows which are not yet frozen. 

To obtain limiting values for the nonequilibrium 
calculations, flow values for the equilibrium case with 
the same pressure distribution have been obtained 
both for real air and for the model gas. The model gas 
values were obtained by the solution of (22), (23), 
(24), and (36) or, actually, with the constant entropy 
condition which is derivable from them and may re- 
place (24). The corresponding values of a are plotted 
in Fig. 4 along with the nonequilibrium values to be 
described. 

Comparing the equilibrium behavior of the model gas 
to that of real air, it is seen that the model gas gives a 
good representation of the equilibrium behavior of 
a (i.e., a) for velocities of 20,000 ft/sec or over for all 
the altitudes considered. In all cases it overestimates 
the change in a@,. At 15,000 ft/sec, the overestimate 
in the change in a, given by the model gas is more ap- 
preciable than at the higher velocities. For engineer- 
ing purposes it may be suggested that when the dif- 


ferences in a, as derived from real air and the model gas 
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are large, the nonequilibrium results to be described 
may be scaled down by the ratio of the equilibrium 
changes of the real gas to those of the model gas to ob 
tain a more realistic view of the nonequilibrium effects 
obtained with the model gas. 

It may be observed from Fig. 3 that the equilibrium 
values of a may be correlated fairly well on the basis 
of a parameter representing the fractional change; 
namely, 


97 


is = as) / Use = Bie) (294 


This trend can be verified further by the use of the 
polytropic exponents for equilibrium air presented by 
Logan and Treanor.*” These yield for constant entropy 


a ww gt (38) 
T~p? 


where y’ and y” are empirical constants given in refer- 
ence 32. Combined with the equation of state, (38) 
yields 


E(2—4) 


] 
lt+aw~p 


(39 
For the present calculations we may take 
vy’ = 1.100 (40) 
vy” = 1.108 
so that 
l-a~ pom (41) 
and 
He — Ae _ (p/pq)?"? — (p:/ pq)? (49) 
Gas =~ Ci 1 — (p laa : 


Eq. (42) is plotted in Fig. 3 and fits the combined 
curves well. 

Now returning to the nonequilibrium calculations, we 
note that various values of L have been chosen for each 
of the altitude-velocity conditions. This fully deter- 
mines the coefficient ¢ of the rate equation (25). The 
resulting values of a and @ are presented in Figs. 4 and 
5. Corresponding typical values of p and T are plotted 
in Figs. 6 and 7. 

It is observed that the values of the characteristic 
length L required to produce noticeable deviations from 
the frozen condition are quite large from a practical stand- 
point. For example, for the 154,000 ft, 25,000 ft/sec 
case, the flow is seen to be substantially frozen for values 
of L < 3 ft. On the other hand, values of Z on the 
order of 200 ft are required to produce conditions close 
to equilibrium. It is indicated that for altitudes 
above 154,000 ft and flight velocities between 15,000 
and 25,000 ft/sec, flows in the inviscid supersonic ex- 
pansion regions of blunt bodies will be essentially chem- 
ically frozen with regard to the molecular dissociation of 
primary components, and may be so treated for theoretical 
analysis and for testing. 

We note that when a complete flow, rather than a 
single point, is considered, there is a degree of arbi- 
trariness in the specification of “‘closeness’’ to either the 
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frozen or equilibrium limiting cases. Here we may 
consider the percent deviation of a or a, from its initial 
value a,, and a flow may be considered frozen if its 
composition changes by less than, say, 5% from the 
initial value. When the equilibrium and _ frozen 
values approach each other as they do at low altitudes 
and velocities, the distinction becomes less meaningful. 
Another parameter considered is the deviation at the 
end of the flow normalized with respect to the maxi- 
mum possible deviation, that is, (@1 — ai) /(Q@ae — ie). 

We recall that L corresponds to the length of a seg- 
ment of a streamline. Its magnitude in the supersonic 
expansion region of a blunt body is on the order of the 
nose radius 7, for a spherical nose, and on the order of a 
fraction of the frontal radius, say r,/4, for a flat-nosed 
body. It would not normally be expected that 
spherical nose radii larger than, say, 5 ft would be 
employed in practice for hypersonic flight above 154,000 
ft. 

For a given Z and velocity the flow becomes more 
frozen at higher altitudes, while for a given L and alti- 
tude it becomes more frozen at lower velocities. These 
trends are due to decreased local densities. More- 
over, for a given altitude and velocity the flow becomes 
more frozen as L is decreased. This is due to the de- 
creased available flow time. 

From Fig. 5 we see that the velocity distributions 
are not altered substantially by deviations from 
equilibrium. In fact, the values of #/tU.,:c do not 
vary appreciably within the range of flight conditions 
under consideration. A check value of the velocity 
obtained for a real air in equilibrium for the 15,000 
ft/sec, 154,000 ft case is shown for comparison at the top 
of Fig. 5. This value is close to that of the model gas 
despite the fact that there is an appreciable difference 
in a, between the model gas and real air in this case 
(see Fig. 4). 

The decreases in temperature (Fig. 7) due to de- 
partures from equilibrium are significant, as may be 
deduced by considering directly the limiting cases of 
frozen and equilibrium flow. The increases in density 
(Fig. 6) stem mainly from the temperature changes; 
the decreases in a make a weaker contribution to the 
density change, as can be seen from the equation of 
state. Since the coefficient of viscosity varies closely 
as 7’? in high temperature air, we may deduce that 
the local Reynolds numbers change roughly in propor- 
tion to 7~*/? due to departures from equilibrium. 
That is, the temperature decreases may result in appre- 
ciable increases in the Reynolds number. In some 
cases the Reynolds numbers increase by as much as 80) 
percent. 

In seeking a general relation between the extent of 
changes in a and values of the dimensionless parameters 
1/y or pa/¢ = pw~CL/u.. in the Prandtl-Meyer flows, 
we have plotted in Figs. 8 and 9 the fractional char ge 
in a, versus each of these parameters. Because only a 
limited number of values of a; have been obtained, the 
faired curves of Figs. 8 and 9 are somewhat hypothetical 
but probably demonstrate the proper trends. 
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Simple correlations are not evident in Figs. 8 and 9, 
although linear variations with In(p.Cl/u.pa) and 
In(p.CL/u..) are shown in the intermediate non- 
equilibrium range. The correlation is not improved 
greatly if p, and u, are used instead of p.. and u., in the 
parameters. Thus an explicit dependence of the a; 
parameter on the initial conditions a, and u,/u. is 
indicated. This question warrants further investi- 
gation. 

Within the limits of the present altitude-velocity 
regime, extreme values of the aforementioned param- 
eters which assure proximity to equilibrium and to 
frozen flow may be estimated, viz: 

as, dh, 


for = (0.95 (frozen): 
Qae = Cie 


< 0.001, ¢ > 1,000 
rt) (43) 


= 0.05 (equil.): 


l 
> 0.01 to 0.1, ¢ < 10 to 100 


r 


¥ 


If we consider a rate parameter analogous to ¢ but 
based on sonic conditions, say ¢,~! = p,’CL/uUgpa, We 
find that for the present flight conditions, ¢/¢, ~ 160 
to 600. Hence the criteria in (43) would be re-expressed 
roughly as follows: for frozen flow, ¢,~' < 0.16 to 
0.6 and yg, > 1.7 to 6; while for equilibrium flow, ¢, 
> 1.6 to 60 and g, < 0.017 to 0.6. These values are 


more reasonable in comparison with unity. 
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Nonequilibrium Flows With Constant Pressure 


A streamline which enters a region of constant pres- 
sure and hence of constant velocity and enthalpy, has 
the opportunity of returning to an equilibrium state 
corresponding to this pressure and enthalpy. The 
extent of the return to equilibrium depends on the 
streamline length L available. 

Several typical calculations have been made for 
flows of this type. Initial states (at s = 0) have been 
assumed corresponding to flows which have emerged 
from the Prandtl-Meyer expansion region with their 
compositions still frozen at the sonic-region conditions. 
Only Eq. (25) must be integrated in these calculations. 
With a new streamwise coordinate defined as 


& = s/Le (44) 


the calculations may be made without prior specifica- 
tion of Ly, whose value may be assigned as a final step. 

The initial conditions for this set of calculations are 
listed in Table 4. Conditions along the streamlines are 
considered only at points less than 100 ft from the initial 
point, since this distance is considered sufficient for 
most practical purposes. 

The results show that, for conditions corresponding 
to altitudes of 200,000 ft and above, and velocities be- 
tween 15,000 and 25,000 ft/sec, the flows are completely 
frozen at stations (values of 5) less than 100 ft from 
the initial point. In the cases of 154,000 ft altitude, 


there are noticeable deviations from the frozen state. 
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TABLE 4 
Initial Conditions for Constant Pressure Calculations 


Alt. 
(ft) 
Uo x 36-9 D p Pd/ Po a tla 71 
15,000 154 2.0559 0.1423 0.872 X 103 0.219 0.296 0.544 
200 2.0922 0.1306 4.630 K 103 0.229 0.291 0.539 
246 2.1362 0.1293 29.762 K 103 0.238 0.275 0.525 
20,000 154 1.4391 0.1326 0.391 0.254 0.504 
200 1.4465 0.1312 0.410 0.238 0.488 
246 1.4523 0.1285 0.429 0.222 0.472 
25,000 154 0.99868 0.1303 0.642 0.223 0.473 
200 0.99939 0.1307 0.669 0.206 0.453 
246 0.99861 0.1323 0.696 0.191 0.437 


Values of a calculated for this altitude are shown in 
Fig. 10. Cross-plots of the values in Fig. 10 are shown 
in Fig. 11. The cross-plots show the fractional change 
in a versus flight velocity at various distances from the 
initial point. Presented for comparison in Figs. 10 
and 11 are the equilibrium values of a which would be 
reached at § > ©. Equilibrium values are shown for 
the model gas and for real air. The air values corre- 
spond to the initial pressures and enthalpies (deter- 
mined by the initial velocities and stagnation enthal- 
pies) given in Table 4. 
At 154,000 ft, the results show that 
Q — As20/ A320 < 0.05, fors < 20ft (45) 
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Although lengths greater than 20 ft have been con- 
sidered in these calculations, it is unlikely that con- 
stant pressure regions extending more than 20 or 30 
ft from the blunt-nose junction will be encountered in 
practice. Indeed, for positions beyond this, and in 
some cases prior to this, we must examine the boundary 
layer flow to discern whether the fluid with stagnation 
conditions corresponding to those behind the body’s 
normal-shock region has not already been swallowed 
by the boundary layer, since streamlines that do not 
pass through the normal shock region have different 
stagnation conditions from those that do. This is a 
Reynolds-number effect which is not treated in the 
present paper. However, we may note that stream- 
lines passing through highly inclined regions of the 
shock will be dissociated to a lesser extent than those 
passing through the normal shock area. Investigation 
of their behavior is left for further work. Only normal- 
shock stagnation conditions are treated in this paper. 

The deviation from frozen flow and the approach to 
equilibrium conditions is most rapid at the higher ve- 
locities, at which the difference in a between frozen and 
equilibrium flow diminishes. We may observe again 
that the model gas overestimates the extent of the 
departure from frozen flow. Therefore the fractional 
change indicated in (45) is probably somewhat larger 
than that which occurs in real air. 

An additional remark may be made concerning the 
conditions which would be found if the pressures along 
the body were decreasing rather than constant. For 
example, such decreases would be found on hemisphere- 
cylinders or blunt bodies with small or negative cone 
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angles in regions beyond the nose. The recombination 
phenomenon dominates over the dissociation in the 
initial region of these flows, as it does in the constant 
pressure flows being described here. This dominance 
persists until the flow approaches the equilibrium 
conditions, at which time the two effects assume equal 
importance. The recombination term in (25) varies 
as p’, and thus is quite sensitive to decreases in pres- 
sure for given initial conditions. Preliminary calcula- 
tions not presented here indicate that even very mild 
favorable pressure gradients lead to flows which are 
appreciably more frozen than comparable constant 
pressure flows. 

Therefore the calculations beyond the Prandtl-Meyer 
region yield the primary result that these flows will be 
essentially frozen, for distances of 20 or 30 ft at least, 
when the pressures are constant or decreasing along the 
body for the flight conditions stated. 


Implications of the Results 


In review, the major sources of empiricism in the 
present analysis are: (1) the estimate of inviscid air 
behavior on the basis of an idealized gas flow; (2) the 
average value of the recombination-rate factor based 
on the experimental results, for oxygen, of Matthews 
and of Byron; (3) the assumed pressure variations; 
(4) the assumption that departures from equilibrium 
are negligible prior to the sonic region of a blunt body; 
(5) the use of stagnation conditions behind strong nor- 


mal shocks. 
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Flow conditions corresponding to altitudes between 
154,000 and 246,000 ft, and velocities between 15,000 
and 25,000 ft/sec, are considered. 

From an analysis predicated upon the aforementioned 
stipulations, we conclude that flows along streamlines 
emerging from the normal-shock region of a hypersonic 
blunt body may be considered chemically frozen with 
respect to molecular dissociation of the primary com- 
ponents in the supersonic ‘‘Prandtl-Meyer’’ expansion 
region, and also beyond it for lengths of 20 or 30 ft in 
constant-pressure or decreasing-pressure zones starting 
from a frozen expansion to horizontal. This implicitly 
restricts the analysis to bodies with small, zero, or 
negative suriace inclinations beyond the nose. 

Although streamlines that penetrate more inclined 
parts of the shock still await analysis, it may be ob- 
served that molecular dissociation phenomena would 
be weaker along these streamlines than along those 
passing through the essentially normal shock region. 
In fact, most of the flow passing through the shock 
region which is appreciably inclined, say, with shock 
inclinations on the order of 45° or less with respect 
to the undisturbed flow, is undissociated. Its thermo- 
dynamics, however, may be complicated by molecular 
vibration. For a thin layer of fluid stemming from 
the region of moderate shock inclination, say, on the 
order of 45°, vibrational processes dominate the thermo- 
dynamics. Vibrational relaxation effects could be 
examined by procedures analogous to those of the 
present paper, but are not treated here. From a 
qualitative standpoint we may distinguish two major 
regions within the shock layer. One will emerge 
from the strong-shock region as dealt with in this paper. 
Assuming translational, rotational and _ vibrational 
equilibrium, the flow will be characterized by values of 
y = (O1n p/O ln p) along streamlines on the order of 1.1 
or 1.2. The other flow regime emerges from the 
weaker-shock region and is characterized by values of 





NOVEMBER, 1960 


SCIENCES- 


y on the order of the undisturbed values, say 1.4. 
Between them is a transitional region previously men- 
tioned. 

Neglecting the effects of vibrational and rotational 
relaxation, we note that the values of y are constant 
along streamlines in otherwise frozen flows, and are 
essentially constant along streamlines (i.e., isentropes) 
in equilibrium flows.*? This fact may be utilized to 
simplify calculation procedures since y may be uniquely 
related to the entropy associated with a given stream- 
line in the shock layer. 

For testing, it may be possible to utilize conventional 
wind tunnels producing flows without the proper dis- 
sociation or without any dissociation, provided that 
the proper values of y are produced in the shock layer, 
and that enough heat is supplied to prevent condensa- 
tion. In this event, force and moment coefficients 
will be adequately simulated. Work is already avail- 
able** ** concerning the simulation of appropriate 
hypersonic air properties, for example, by the use of 
gases other than air and particularly with products of 
chemical reaction. 

The equilibrium values of a,, calculated in the 
Prandtl-Meyer region and correlated in Fig. 3, em- 
phasize the fact that the constant polytropic exponents 
vy’ and y” of reference 32 may be used along the stream- 
lines (lines of constant entropy in equilibrium or frozen 
flow) not only to relate the pressure, density and tem- 
perature, but also to yield the atom concentration 
fraction a,. In fact it may be suggested that properties 
of equilibrium flow fields may be obtained approxi- 
mately by first assuming constant-composition flow 
to calculate pressure distributions and stagnation condi- 
tions along the streamlines. Then these pressure dis- 
tributions can be employed together with the constant 
polytropic exponents to give approximate equilibrium 
values of temperature, density, atom concentration, 
velocity and electron density, and corrected stream- 
line shapes. Of course, this procedure is only a simpli- 
fied version of the streamline calculation procedure 
employed previously in this paper for the nonequilib- 
rium and equilibrium calculations. Its simplification 
rests in the use of the polytropic exponent correlation 
of reference 32. 

The streamline analysis procedure used in this paper 
may also have utility for experiments. In an experi- 
ment, a converging-diverging nozzle would be designed 
to produce the desired pressure distribution along a 
streamline of desired length. If it is required to in- 
vestigate real gas effects, such as the extent of non- 
equilibrium dissociation along the streamline, it is 
necessary to duplicate true stagnation conditions by 
means of either shock tubes or plasmas. 

In further work along the lines started in this paper, 
it would be interesting to investigate other streamline 
conditions, such as the flows emerging from inclined 
shocks. Also, complete flow fields could be examined 
for nonequilibrium effects. It would also be useful to 
utilize more exact representations of the chemical 
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reaction phenomena and to compare the results with 
those presented here. 

Finally, we note that the flows considered here are 
already completely frozen at altitudes for which rare- 
fied gas phenomena of a fluid-mechanical rather than 
fluid- 
effects include increases in 


a chemical nature come into play. These 
mechanical rarefaction 
vorticity in the shock layer, boundary-layer thickening, 
and shock thickening, and begin to manifest them- 


selves at altitudes roughtly above 150,000 ft. 


Conclusions 


One-dimensional, inviscid nonequilibrium flows of a 
two-component model gas have been studied for 
prescribed pressure variations and an average reaction 
rate based on recent data for oxygen recombination. 
These flows have been interpreted in relation to the 
flow along streamlines around blunt hypersonic bodies. 
Assuming equilibrium conditions in the subsonic region, 
it is estimated that the flow in the initial supersonic 
expansion region, which is approximately of Prandtl- 
Meyer character, will be chemically frozen with respect 
to the molecular dissociation of the primary com- 
ponents under the hypersonic, high-altitude flight 
conditions considered. The flight conditions consist 
of flight velocities between 15,000 and 25,000 ft/sec at 
altitudes between 154,000 and 246,000 ft. Further- 
more, on bodies of small surface inclination beyond the 
nose, the flow will continue to be effectively frozen 
for at least 20 ft downstream of the nose. These con- 
clusions may lead to the simplification of procedures 
for theoretical calculation and testing. 

The problem of distinguishing a dimensionless 
length-reaction rate parameter characterizing the ex- 
tent of departures from equilibrium or from frozen 
behavior in the flow fields of interest here has been 
discussed. For the supersonic expansion (‘‘Prandtl- 
Meyer’’) region a length-reaction rate parameter 
Ga! = pa*CL/ugp, is defined on the basis of sonic 
conditions as an alternative to the analogous param- 
eter ¢~! based on flight velocity and density. On this 
basis, the following requirements in the supersonic 


region may be stated: 


for frozen flow: Ga? < 0.16 to 0.6 
for equilibrium: ¢,~!> 1.6 to 60 
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An Optical Boundary-Layer Probe’ 


JOHN W. DAIBER* 


Cornell University 


Summary 


The applicability of the schlieren-photomultiplier technique to 
obtain quantitative density measurements in the laminar bound- 
ary layer induced by a traveling shock wave in a shock tube is 
investigated. Tests were conducted at a Mach number of 1.58 
so that the data could be compared with the exact theoretical 
solution tabulated by Mirels. The data obtained are in good 
agreement with the theory if the distance of the light beam above 
the floor of the shock tube is adjusted to fit the theoretical curve; 
this would not be necessary if a larger shock tube were used. 
Values of the transition Reynolds number were also determined 
which are slightly less than those found by Martin using an inter- 
ferometer. It is shown that this technique is sensitive enough 
to detect changes in density that are only 0.006 per cent of at- 
mospheric density. 


Symbols 

A = signal strength from the light screen 
I) = light intensity predicted by geometrical optics 
I; = light intensity predicted by physical optics 
n = index of refraction 
t = time 
T = absolute temperature 
Us = shock velocity in lab coordinates 
Uy = fluid velocity in lab coordinates 
V = voltage 
x,y = steady, shock-fixed coordinates, y normal to plate, x 

positive ahead of the test section 
Ye = position of upper edge of light beam 
y = position of lower edge of light beam 
ri) = boundary-layer thickness 
€ = angular change in beam direction 
@ = acute angle between shock and light beam 
“a =6= Viscosity 
p = density 
p) = reference density at one atmosphere and (°C 
( ); = condition before the shock 
( )2 = condition after the shock 

Introduction 


HEN A SHOCK WAVE advances into a stationary 
fluid in a shock tube, it sets the fluid into uniform 
motion at a constant temperature and pressure. Because 
the fluid has viscosity its motion will create a boundary 
layer on the walls of the tube. This boundary layer 
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has an effective leading edge which is always at the 
shock. This unsteady problem can be reduced to a 
steady problem by transforming to a coordinate system 
in which the shock is stationary. Several authors? * 
have made the above transformation and other simpli- 
fying assumptions to obtain theoretical solutions to the 
laminar boundary-layer problem. 

It is the purpose of the present study to develop an 
experimental method for obtaining accurate, quantita- 
tive density data within the boundary layer. A 
method is investigated which combines the sensitivity 
of the schlieren with a photomultiplier tube to obtain 
absolute values of the local density. Some experi- 
mental data are obtained, and the method of interpret- 
ing these data is presented. The reduced data are then 
compared with the theoretical solution tabulated by 
Mirels, and the validity of the theoretical assumptions 
is verified. 


Description of Apparatus 


The shock tube was constructed from 13-gage, 
welded, cold-rolled, rectangular steel tubing, with the 
internal burr removed. The internal tube dimensions 
were 1.31 and 0.81 in. normal and parallel to the light- 
beam direction, respectively. The driven section was 
14 ft. long with 1/2-in. thick circular plate glass test 
section windows 11'/, ft. from the diaphragm. The 
driver section was 16 in. long with a piloted flange that 
mated with a 1/8-in. recessed flange on the driven sec- 
tion. An aluminum foil diaphragm was placed between 
the flanges which were then held in place by two lock- 
jaw wrenches. A leak proof seal was obtained at the 
diaphragm by having a neoprene gasket on the driven 
flange and a 0.003-in. protruding rib on the driver 
flange. 

In the shock tube was a flat plate of plate glass ex- 
tending 28 in. ahead of the test section which reduced 
the test section height to 0.96 in. 

A shield before the test section divided the parallel 
light beam from the source, an Osram HBO 107/1 
lamp, into three beams. After the first two beams 
passed through the test section, stops were placed in 
them to make use of the light screen method of measur- 
ing shock speed. These two beams then entered a 
photomultiplier tube whose output was viewed on an 
oscilloscope by way of a 120 cycle filter to remove modu- 
lation of source intensity. This scope was triggered 
internally on the signal from the light screen. 

The third beam was 1/8 in. high and was focused on 
an horizontal knife edge. This beam then entered a 
photomultiplier whose output went directly to a second 
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scope which was triggered by the delay gate from the 
first scope. 


Theory 


(a) Flow Behind Shock 


Consider a flat plate of zero thickness to be placed in 
the shock tube at zero angle of attack. When the 
shock wave passes onto the plate a boundary layer will 
build up on the plate. If the shock wave is traveling 
with velocity U, and the fluid after the shock has 
velocity ly, then particles initially located above the 
leading edge of the plate have in a time ¢ moved over 
the plate a distance Uot, Fig. 1. The particles in 
region A originated ahead of the plate’s leading edge. 
Particles which were initially above the plate are in 
region B and are affected only by the shock.’ 

The relative size of the interaction region A decreases 
with decreasing Mach number. At 28 in. from the 
leading edge of a flat plate, the time difference between 
arrival of an J = 1.58 shock and the arrival time of 
the interaction region is more than ten times that re- 
quired for transition from laminar to turbulent flow in 
the boundary layer when the pressure before the shock 
is 1'/. in. Hg (see reference 4). Thus, measurements 
may be made on the plate without consideration of plate 
leading-edge effects. 


(b) Density Measurement 


With a schlieren system, the change in illumination, J, 
at an image point is directly proportional to the product 
of density gradient and illumination at the correspond- 
ing point in the test section, when the flow is two- 
dimensional.® 


I « (I,/I)(dp/dy) (1) 


The conditions of the experiment were chosen so that 
the deflection of the light beam through the boundary 
layer constituted a negligible displacement of the beam 
compared with the boundary-layer thickness. 

By placing a photomultiplier tube directly behind the 
focal point of the schlieren system, quantitative density 
measurements can be obtained.® This is possible be- 
cause the illumination from the schlieren system is pro- 
portional to the density derivative and the current 
produced by a photo-multiplier is proportional to the 
total illumination or the integral of the density deriva- 


tive. Hence, the voltage above the base value 6V 
produced across the photomultiplier’s load resistor is 
given by 


= -I,dp 
6bV « —— dy (2 
J, Indy ~ , 


As mentioned above, a shield was used to define a 
beam 1/8 in. high through the test section. Actually 
the shadow of this shield is not in sharp outline through 
the test section, but consists of a system of dark and 
light bands in the immediate neighborhood of the edge. 
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These bands constitute the well-known Fresnel diffrac- 
tion pattern. 

The upper edge of the beam is in the free stream where 
blurring is of no consequence because the density 
gradient is zero there. Therefore the upper limit can 
be replaced by y,._ If there were no diffraction effects, 
the intensity would be a step function at yo, in which 
case Eq. (2) integrates to 


i Pu — Pue (3) 


If a constant density gradient is assumed for the 
boundary layer, then, because interference effects are 
conservative of light energy, Eq. (2) will again reduce 
to Eq. (3). 

Once the constant A has been experimentally deter- 
mined for a particular case, as, for example, by the 
wedge technique as suggested in reference 6, and carried 
out in references 7 and 8, a measurement of the photo- 
multiplier voltage gives the local density relative to the 
free-stream density. 


(c) Density Profiles 


By assuming similar profiles, a relation is obtained 


of the form 
P — Py, [y 
=f 4) 
Po : (7) \ 


The boundary-layer thickness 6 depends on the ex- 
ternal steady flow velocity, U2, the kinematic viscosity 
evaluated at the wall, v,, and the distance from the 
shock, x. The form of this relation is taken from 
steady-state flow, where 


6 « V 2x4 U2 (5) 


These variables are not directly measurable, but can 
be reduced to quantities which are. Because the pres- 
sure across the boundary layer is constant, the equation 
of state for a perfect gas can be used to give 


Vo = Mel o/ Pel 

The subscript w means that the quantities are to be 
evaluated at the wall after the shock passes. However, 
it is a very good approximation to assume that the wall 
temperature after the shock is the same as before the 
shock. This is valid for small testing times, since the 
laminar boundary layer has a low heat-transfer coef- 
ficient and the wall has a large heat capacity.® 

The steady flow velocity after the shock is related 
to the velocity before the shock by the conservation of 
mass. 


U2/x = (pi/pe)(Us/x) 


Laboratory time is determined by the speed of the 
shock, so 


Us/x = pi/ pe 


Thus, Eq. (4) can be expressed as 
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=] Vo aie (6 
Po , ( \ 2uil | ” 


In the present method of obtaining experimental data, 
the location of measurement above the wall is fixed at 
yo. The function is obtained in one run from the time 


variation. 

In each run the shock speed is measured; then, by 
knowing the initial temperature, the Mach number can 
be calculated. The Rankine-Hugoniot shock relations 
can then be used to find the density and temperature 
ratio across the shock. The initial shock tube density 
is found by measuring the pressure and the temperature 
prior to breaking the diaphragm. 


(d) Fresnel Diffraction 

The assumption of a constant density gradient in the 
laminar boundary layer is very poor; therefore Eq. 
(2) was graphically integrated for the beam at a fixed 
location of 0.004 in. above the plate at five different 
A typical curve of the integrand in Eq. 
The integral was then 


testing times. 
(2) vs. y is shown in Fig. 2. 
evaluated using a planimeter for both the actual varia- 
tion and the ideal variation when the intensity is a step 
function. Table 1 summarizes the results. 

For early times the percentage error is large. 
ever the actual increment difference could only be de- 
tected with an instrument of low noise level. 

Thus, for the investigation of a laminar boundary 
layer the Fresnel diffraction is most definitely present, 
but its effect on the data is completely negligible (for 
shock Mach numbers of 1.58) if the distance yo is at 
least 0.004 in. 


How- 


Measurement of Shock Speed 


Through a shock wave the density increases so the 
index of refraction after the shock is greater than the 
index of refraction before the shock. When a narrow 
light beam is incident with a grazing angle @ on an 
interface where the index of refraction increases, the 
light beam is refracted, an angle ¢ relative to its initial 
direction. If a shield is placed in the initial beam aft 
of the location where the shock will occur, then as the 
shock crosses the beam, light will be deflected off the 
shield. By placing a photomultiplier tube behind the 
shield, the deflected light will cause a rapid increase in 
signal. This then gives a simple method, called the 
light screen technique, of measuring the shock speed by 
using two such beams and recording the time between 
deflections. 

The angle of deflection is quite small, and the initial 
and deflected beams are still partially coincident for 
the shield a reasonable distance from the tube section. 
The strength of the signal for a fixed location of the 
shield depends on the product of the amount of light 
refracted, given by Fresnel’s equation,'® and the angle 
of refraction, given by Snell’s equation.’ The signal 
strength A is thus given by 








SCIENCES—NOVEMBER, 1960 
TABLE | 
Predicted Photomultiplier Voltage With and Without Diffraction 
Effects, JJ = 1.58, vw = 0.004 in., P; = 1.50 in. Hg 
t (microsec. ) (6V/ VK )actual (6V/VK )idea 
BS 0.00086 0. 00027 
8 0.00234 0. 00206 
5 0.00581 0.00580 
30 0.00944 0.00968 
85 0.01741 0.01728 


m\" . . cos [arecos (71,72) cos 6 — 6] . 
Aa sin? 20 —— (7) 
nN» sin” [arccos (7/12) cos 80 + 6] 


This equation is plotted in Fig. 3. As can be seen, 
there is a definite maximum to the signal strength, the 
location of which depends on the relative index of re- 
fraction for the shock. 

In the present study the relative index is the order of 
1.00005, so the optimum grazing angle is 30 min. Con- 
sequently, the light source was set so the beam would 
make 1/2° angle with the normal to the shock tube 
walls. The two shields are placed in beams one and two 
at the maximum distance possible from the test section. 


Error 


The space resolution due to a finite source size and 
bending of the light beam by the density gradient was 
less than 0.0005 in. This gave a percentage error of 
13 per cent with a yo of 0.004in. This percentage error 
could be greatly reduced by testing at lower pressures 
so that the boundary layer would be thicker, allowing a 
greater y) to be used. The largest error in locating yo, 
which was 100 per cent in the present experiment, is 
caused by the light beam not being parallel to the flat 
plate. This could be corrected in a shock tube with a 
test section of greater height by using the combination 
of a first surfaced glass plate and a penta-prism as an 


autocollimator. 


Results 


The success of the light-screen technique in measuring 
Mach number can be seen in Fig. 4. The two signals 
have a rapid rise time which makes it possible to meas- 
ure the distance between them accurately. 

The data presented in Fig. 5 was reduced by the 
method outlined earlier and is shown in Fig. 6. The 
form of the continuous experimental curve closely 
resembles that of the theoretical curve, but there is a 
large displacement of one curve with respect to the 
other. This discrepancy is due to the error in deter- 
mining the value of yp. The true value of yo can be de- 
termined by adjusting one point to the theoretical 
curve. By doing this the line of experimental data 
agrees closely with the theoretical curve. 

As was shown in the theory, the density has a square 
root dependency on the time. When the experimental 
traces from the oscilloscope are plotted on log-log paper 
the curve should be linear up to the time that transition 
from laminar to turbulent flow takes place within the 
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boundary layer. So, by measuring the time at which 
the trace broke from a square root dependence on time, 
the schlieren-photomultiplier was used to record the 
transition distances. The Reynolds number at transi- 
tion based on the wall temperature which were obtained 
are in close agreement with those obtained by Martin,‘ 
as can be seen in Table 2. The estimated error in the 
transition Reynolds number is at most 10 per cent and 
is determined only by the noise level. 


Conclusions 


The applicability of the schlieren-photomultiplier 
combination to obtaining the transition Reynolds num- 
ber and the value of local density within the boundary 
layer has been demonstrated. Density differences as 
small as 1/2 per cent of atmospheric density were easily 
detected. As shown in reference 8, by using a d.c. source 
with no modulation of the light intensity, density 
changes amounting to 0.006 per cent of atmospheric 
density could be detected. This would then allow a 
very careful and accurate study to be made of the 
boundary layer by a method which probes the flow 
field optically. 


References 


1 Hollyer, Robert N., Jr., A Study of Attenuation in the Shock 
Tube, Eng. Res. Inst., University of Michigan, July, 1953. 

2 Mirels, Harold, Laminar Boundary Layer Behind Shock Ad- 
vancing into Stationary Fluid, NACA TN 3401, March, 1955. 


Plastic Bending (Continued from page 820) 


SCIENCES—NOVEMBER, 1960 





TABLE 2 
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Closed-Loop Flip-Flop Control Systems 


MILTON ROGERS* ann GEORGE SHAPIRO** 
Air Research and Development Command, USAF, and Westinghouse Electric Corporation 


Summary 


On-off, or flip-flop control systems, as they are called, are 
assuming an ever-increasing, important role in the development 
of advanced missiles and rockets. The primary limitation to the 
inmmediate acceptance and use of flip-flop control systems by 
the engineering community lies in the fact that such systems are 
essentially nonlinear. The development of such systems, there- 
fore, departs from the body of knowledge and experience gained 
by the engineer in designing linear control systems, which have 
been so successful to date in advancing the horizons of controlled 
flight and automation. 

Here the basic concepts underlying the synthesis and analysis 
of closed-loop, flip-flop control systems of the electromechanical 
type are developed. Generalized design and analysis charts are 
presented to aid the engineer in gaining quickly an understanding 
of the relative importance of various elements in the system to the 
system’s performance and stability. The importance of con- 
sidering not only control-system dynamics, but also the damping 
provided by the airframe in developing practical flip-flop control 
systems is forcefully demonstrated to the engineer and scientist. 
Analog computer data, obtained as part of the present study, are 
used to substantiate analytic approaches and to advance the 


understanding of such systems. 


Symbols 


For convenience in design, the solution is presented in terms 
of the following quantities which can be determined readily in the 
laboratory or from component manufacturers’ specifications. 


te = total time for maximum control torque in one direction 
to become maximum control torque in the other 


direction, sec. 


tp = presumed delay time in the system before the actuator 
or system responds after an energizing signal, sec. 

fy = control dwell time in a maximum control torque condi- 
tion, sec. 

t, = lead time, defined as —6/6 at signal trip, sec. 

6 = roll angle, rad. 

K = aerodynamic damping in roll = —C),(b?/4)pSV, ft. Ib. 
sec. 

Ci, = aerodynamic damping coefficient in roll = 0C,/0(6b/2V) 

q = dynamic pressure = 1/2 pV?®, lb./ft.? 

p = ambient air density at the altitude under consideration, 
slugs /ft.8 

V = forward velocity, ft./sec. 

S = characteristic area, ft.? 
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b = characteristic span, ft 

I = moment of inertia of the vehicle in roll, slugs ft.? 
L = applied rolling torque, ft.lbs 

k = normalizing time parameter K|/I, sec.— 

T = normalized time = kt 

x = normalized roll amplitude 0K2/IL, rad 


Dot and prime indicate differentiation with respect to ¢ and rf, 


respectively 


Introduction 


| apie CONTROLS are usually of the proportional 
type, wherein the applied correction moment is 
proportional to an error in the ‘‘attitude,” or they may 
be of the flip-flop (‘‘on-off’’) type in which a maximum 
corrective force is applied whenever an error of sufficient 
magnitude is sensed by some measuring device. 

Proportional systems, using reaction controls, have 
the dubious advantage over on-off systems in that linear 
mathematical techniques may be used for both analysis 
and synthesis. On-off systems have advantages over 
the more common proportional systems in that they 
weigh less, are simpler to build and more reliable 
in use since they require relatively few components, 
and need less power to operate. They are, there- 
fore, particularly attractive as airborne controls, wherein 
a limited source of power is available to operate the 
system and where low weight and high reliability are 
demanded. However, on-off type controls are non- 
linear, and their successful incorporation into aero- 
nautical technology will require a major departure from 
analysis techniques used so successfully by aeronautical 
engineers in the past. 

Optimal control systems are those which provide 
simultaneously for minimum error (either static or 
dynamic), minimum response time, stability (including 
reliability) of operation, minimum power consumption, 
and minimum complexity. Any practical design is a 
compromise between these conflicting considerations 
and the relative importance given to these factors by 
the control system designer will be determined by the 
performance requirements of the entire system. The 
ability of the designer to properly synthesize a control 
system, however, depends on how much is understood 
of the relationships between the parameters of a system 
and their effects on system performance. 
Thus, proper synthesis is contingent upon adequate 
analysis. 

It is the purpose of this paper to analyze in a realistic 
fashion the performance of a general class of airborne, 
closed-loop, on-off control systems. The contributions 
of various control and airframe characteristics including 
aerodynamic damping to the overall system perform- 


overall 
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Fic. 1. Block diagram of closed-loop, on-off roll control system 


ance are delineated and discussed. Analytic procedures 
for engineering use are derived. Approximations and 
design charts are presented to aid the engineer in gaining 
quickly an understanding of the system performance. 
Since optimal control systems are dependent upon 
specific performance requirements, optimization of a 
system is beyond the scope of the generalized study of 


this paper. 


Description of System 


Although an aerial vehicle is free to respond to a 
control input in pitch, yaw, and roll, only the rolling 
motion is considered in this analysis; that is, the on-off 
system analyzed is considered to be a single-axis con- 
trol. This simplification of the intricate overall prob- 
lem of vehicle dynamics does not violate any physical 
reality if the rolling amplitude is small, and it permits a 
straightforward evaluation of the contribution of the 
various system parameters to the system performance. 
The analytical procedures developed are equally 
applicable to the other motions, such as pitch or yaw. 

It is assumed that the aircraft is subjected to a rolling 
moment that is directly proportional to the instantane- 
ous roll-control displacement or torque, and that the 
aircraft reacts to this roll control moment in a manner 
dictated not only by inertial considerations but also by 
the aerodynamic damping in roll provided by the air- 
frame. The roll damping, as used, is defined as a 
moment opposing the rolling motion and is considered 
to be proportional to the angular rolling velocity. In 
mathematical terms the aircraft, or missile, considered 
is assumed to behave under the influence of the control 
as a linear, second-order system with constant co- 
efficients. 

It is well known that the aerodynamic damping of a 
more conventionally stabilized and controlled vehicle 
exerts a predominating influence on the vehicle sta- 
bility, and also on vehicle response, rolling amplitude, 
and rolling velocity. This damping depends, of course, 
on the flight speed and altitude, as well as on the ex- 
ternal shape of the airframe, including tails, body, and 
wings (and especially on the size of these aerodynamic 
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surfaces). Neglect of damping in the analysis and 
synthesis represents a major simplification of the 
mathematics involved! and may be physically proper 
in certain cases. It is one of the purposes of this paper 
to delineate the importance of considering this aero- 
dynamic damping in the design of on-off control sys- 
tems and to indicate those design conditions in which 
it is ‘“‘safe’’ to neglect the damping. Thus it is hoped 
to put the synthesis of airborne on-off control systems 
onto a more rational physical basis than would appear 
to be the present case. 

In any physically proper mathematical model of a 
system such as a roll control, it cannot be assumed that 
maximum torque, in either direction, can be obtained in 
zero time. The rolling torque, be it from an aileron or a 
reaction control, takes a finite time to build up to 
maximum effectiveness in one direction, after having 
been at maximum effectiveness in the other direction. 
There are several ways in which this torque build-up 
time can be considered. One method could be to 
assume that the torque build-up acts as a time delay in 
the system.' This would represent a simplification in 
that only a step pulse need be considered. However, 
the control is providing a torque, which gives rise to an 
acceleration, during the build-up and delay period. 
Thus, although this model perhaps provides a mathe- 
matical simplification, it does not represent the physical 
reality, and its use in any general study may lead to a 
loss of physical significance. 

For the analysis considered herein, therefore, the roll 
control is not assumed to arrive at the condition of 
maximum torque in zero time, but is assumed to take a 
finite time to build up to, and decay from, this maxi- 
mum effectiveness. The effectiveness of the control 
during the build-up and decay is considered to increase 
linearly with time so that the control torque is mathe- 
matically represented by a straight line ramp function 
followed by a flat or plateau, rather than as a pure step 
function. The assumption that constant control 
torque build-up or decay rates are attained instantane- 
ously is not perfect. It is relatively easy to formulate 
empirical equations which will closely approximate any 
given torque build-up or decay. However, these 
equations become unwieldy to use in a generalized 
solution and their use could not contribute any more to 
the general understanding of the relationships between 
the various parameters than the assumed straight-line 
ramp function. 

Any mechanical or electromechanical system cannot 
respond to an error signal instantaneously, that is, 
without some time delay in the response to a control 
error input. Physical control systems are composed of 
valves and other moving components, all of which have 
inertia and friction and therefore take a finite time to 
respond to a signal to move. The system may also 
contain linkages and gear trains, and backlash may thus 
be involved in the response. This backlash also may be 
considered as an effective time delay in the system. 
For the purpose of this analysis, all delays in responding 
to a signal are lumped together as a total system delay 
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time, fp, with no attempt made to delineate which 
particular component gives rise to which portion of the 
total time delay. 

A block diagram of the system considered is given in 
Fig. 1. Unlike proportional control systems, a practical 
on-off control system must be of the closed-loop type, 
whether it be by position feedback only, or by velocity 
feedback, or by any combination of the position error 
and its derivatives. As can be seen from the block 
diagram, both position and velocity are used in the 
feedback circuit. Following conventional engineering 
practice, the position and velocity signals are combined 
linearly. 

Jt is to be noted that a properly operating flip-flop 
control system is, by definition, a neutrally stable sys- 
tem whose steady-state operation is characterized by a 
limit cycle. The steady-state control response of flip- 
flop control systems is shown as a time trace of the 
control torque in Fig. 2. The various system time 
constants are indicated for clarity. 

Neither hysteresis nor a dead zone are considered in 
the control sequence analyzed. These factors can be 
included in the analysis without radically changing the 
fundamental approach indicated. An analysis of their 
effects on the system performance increases the tedium 
and detail of the analysis and does not appear to in- 
crease the fundamental understanding that can be 
extracted from the analysis given. 

By suitable determination of the various system 
parameters chosen, the analysis, as given, is equally 
applicable to ailerons, elevons, or on-off reaction-type 


jets. 


Analytical Methods 


Mathematical Model 

In the previous section, the vehicle under consider- 
ation was described as behaving as “‘a linear, second- 
order system with constant coefficients.” Thus, 


T6(t) + Ké(t) = Lit) (1) 


Reference to Fig. 2 shows that, for time measured in 
each case from the break points of the control torque 


cycle, 
Lit) = —L he aig 
= —L+ (2Li/tz) OStS te ie 
= +L “i tle is 
= +L — (2Lt tr) 0 < t < tr 


and at the trip points (¢ = fy in the first and third seg- 


ments of Eq. (2) above): 
6+ t,6=0 


It would appear that specification of the aircraft 
characteristics requires the use of two parameters such 
as K/I and L/J. However, by the proper normaliza- 
tion of the units of time and angle these parameters may 
be eliminated. Put 
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Fic. 2. Time trace of control torque 


T K\t/1 kt 
x K*arh Rl6 L 
D,(t) Li(t)/L 


Then (1) and (2) become: 


” K ; 
+ ange D,{|t(r) D(r) 
K 
D(r) — |] OS 7rS tH + Tp 

—1 + (27/Tp) O<4r< te 

+ Os +r S ta + tp 

+] — (27/TRpR) O< rT < te 
where 

tH = kty, to Rtp, Tr Rte, tr = kt, 
K|\/K = +11f K> 0, A\/K —-1ifA<0O 
and at the trip point (r = ry in the first and third seg- 
ments): 
x + rx’ 0 


Thus the aerodynamic characteristics are removed by 
the normalization process. Hereafter, except in the 
numerical example discussed later, only the normalized 
form of the system will be considered. 


Mathematical Techniques 

Various techniques for the analysis of on-off systems 
have been used. These are summarized in reference 2 
which also contains an extensive annotated bibliog- 
raphy. Generally, these techniques may be character- 
ized as frequency-domain methods, phase-plane meth- 
ods, or time-domain methods. 

The describing function technique is a frequency- 
domain method. This is an attempt to extend to non- 
linear systems the methods applicable to linear systems. 
For a periodic input function, the output of the non- 
linear element may, in general, be expanded in a Fourier 
series. The remaining linear elements of the system 
will usually attenuate the higher harmonics of the out- 
put of the nonlinear element. Therefore, it is assumed 
that this element is adequately described by the funda- 
mental component of its output. As is pointed out in 
reference 2, in a number of cases this has led to errone- 
ous results. 

The phase-plane trajectory of the solution is a plot of 
x’ vs. x. This approach is particularly well suited to 
second-order equations in which time does not appear 
explicitly. When, as is true for the system under 
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Fic 3. Schematic of analog simulation. 
consideration, several functional forms are required for 
an analytic description of the system, matching of the 
end points of segments of the phase-plane trajectory is 
required. For a particular system, this may often be 
performed graphically; however, it is difficult to obtain 
general analytic results. Time parameters, such as our 
Tr and rp, particularly complicate such an analytic 
treatment. The major utility of phase-plane methods 
lies in the understanding obtained of the transient 
response of the system. 

Time domain approaches are usually possible only 
when explicit solutions may be obtained. This is 
always true for a piecewise linear system such as is being 
considered herein. These approaches also require the 
matching of end point conditions. An elegant, al- 
though highly formal, method of simplifying this pro- 
cedure is due to Kahn and is also described in reference 
2. However, Kahn’s method does not appear to be 
applicable to the case in which the forcing function, 
D(r), is not a step function. In the present case, D(r) 
includes ramp segments. 

In Appendices A and B an explicit time-domain 
approach, supplemented by occasional recourse to the 
phase plane, is used. These methods, though not at all 
sophisticated, lead to closed-form relations between the 
system parameters, lead time, delay time, rise time and 
the limit-cycle parameters, peak amplitude and dwell 
time. (The limit-cycle period is 2[rp + tr + Ty.) 
These relations are necessarily satisfied for any limit 
cycle, but have not been shown to be sufficient to ensure 
the existence of such a limit cycle. Thus the major 
value of this technique, as applied here, is the quantita- 
tive description of the steady-state behavior of the 
system. 
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Analog Computers 


A powerful method available to the control-system 
engineer is the use of computers to aid in the solution of 
intricate systems of differential equations. The most 
useful of the computers for dynamical analysis and 
control system synthesis is the analog computer. 

The Westinghouse Air Analog Computer 
Facility was used to verify the analysis made. A 


Arm 


schematic diagram of the analog simulation is shown 
in Fig. 3. Close examination of this schematic will 
reveal several facets of the simulation problem that 
deserve special mention. 

As with any engineering solution to a problem, the 
final ; roduct is a compromise between several conflict- 
ing factors. The problem of analog simulation of a 
generalized system is in some respects easier, but in 
many aspects more difficult, than the problem of simu- 
lation of a particular control system whose components 
can be used as part of the analog. 

In the interest of simplicity and certainty of perform- 
ance, it was decided to avoid the use of relays, realizing 
that not all of the data obtained would be useful in 
checking the analytical results. For example, in some 
cases the circuitry used to provide the desired delay 
and rise times could not have sufficient time to reset 
could become 


between cycles. In other cases, ty 


negative and the simulated system could “‘chatter.” 

During analog runs the operation of the circuits pro- 
viding the rise and delay time was monitored, so that 
improper simulation could be observed. 


Results 


The analytical work followed three different yet 
interrelated paths: a theoretical development and 
analysis, analog studies to verify the theory and obtain 
additional data, and a cursory engineering study to 
demonstrate the practical significance of the results. 
For convenience and emphasis the results of the study 
are summarized and correlated below, while the de- 
tailed mathematical derivations of the final expressions 
for dwell time and limit cycle amplitude are given in 
Appendices A and B respectively. 


Dwell Time 

The cyclic frequency of the system during the limit 
cycle is determined by the sum of three time terms—the 
rise time, delay time, and dwell time. The relationship 
between these parameters and the system lead time is 
given by the transcendental Eq. (Al1) of Appendix A. 
This relation gives r, explicitly as a function of the 
normalized parameters ty, Tp and rr, whereas ry is 
obtainable only numerically or graphically. To aid the 
engineer in obtaining rapid estimates, the relation is 
presented graphically in Figs. 4 and 5. 

It may be seen from Fig. 4 that the effect of rise time 
is only significant for small dwell time. For large re 
and ry, the expression 
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y= = = 


is approximately valid, although for ty < 2 it under- 
estimates the value of rz. As expected, increasing the 
lead time 7, decreases the dwell time of the control. 

At tr = 1.6 the relationship among the parameters 
changes markedly, both qualitatively and quantita- 
tively. Figure 5 shows the relation between 7, and ry 
for several values of rp when rz = 0.5. For the larger 
values of dwell time the relation is similar to that 
obtained for large rise time, but as lead time is increased 
the dwell time approaches a generally positive limit 
which depends on the delay and rise times. In order to 
reduce dwell time below this value, negative r,, that is, 
lag rather than lead, is required. This effect cannot 


occur for rr > 1.6. 
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Fic 6. Stability boundary 

Since the upper branches of the curves of Fig. 5 
eventually (as ry — ©) become straight lines of slope 
—0.5, it follows that for certain small values of rise time 
there are values of rp and r, which may theoretically 
result in two different limit cycles. As noted in Appen- 
dices A and B the relationships discussed above are 
necessary for the existence of limit cycles but have not 
been shown to be sufficient. Assuming both limit 
cycles to exist mathematically, the question of their 
physical existence, that is, their stability, may also be 
raised. Fliigge-Lotz* in considering a similar problem 
also found the possibility of two limit cycles. She 
states that it can be shown that only the shorter limit 
cycle is stable. It appears from our work that only the 
longer limit cycle is stable; this dissimilarity may stem 
from the difference in the problem considered. 

Figure 5 also shows all of the pertinent analog- 
computer data. Obviously, the agreement is excellent. 
The only visible discrepancies occur at rp = 0; even 
there the agreement is good, especially in view of the 
inherent small delay in any analog simulation. Un- 
fortunately, no analog runs were made with input con- 
ditions which theoretically might result in double limit 
cycles. Thus, the analog results shed no light on the 
question of the coexistence of two limit cycles. A 
phase-plane plot for the case tr, = —0.25, rp = 1.0, 
Tr = O indicated the existence of an unstable inner 
(small-amplitude) limit cycle, as well as a stable outer 
limit cycle. However, the inherent inaccuracies of the 
graphical technique preclude any firm conclusion. 

Figure 6 obtained from Eq. (A12) shows the system 
“stability” boundaries defined as ry = 0. In some 
cases these curves have two branches when rz is suffi- 
ciently small. When two branches exist, one ‘‘stable”’ 
limit cycle is possible in the region below the upper 
branch, while in the region below both branches of the 
curve two “‘stable’’ limit cycles are possible. 

The term ‘‘stable’’ is used in an unconventional 
manner in the preceding paragraph. The intent is to 
describe a system which operates ‘‘as expected,”’ that is, 
the trip signals, x + r,x’ = 0, occur sufficiently far 
apart in time that after one signal the control can com- 
plete its motion to the alternate ‘‘hard over’ position 
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(after allowing for system delay time) before the next 
trip signal. The ‘“‘unstable’’ limit cycle is analogous to 
the chattering of a relay. In control systems, the 
physical behavior depends on the mechanical detail of 
the equipment; in some cases (reaction controls) the 
control may actually ‘‘chatter,’’ while in other systems, 
the control must travel to the “hard over’’ possition and 
can return only after the requisite delay. In such 
systems, the control actually operates on the “‘stability”’ 
boundary with an effective 7, given by Fig. 6 rather 
than by the design value. 


Limit-Cycle Amplitude 

The amplitude of the limit-cycle oscillations is shown 
in Figs. 7 and 8 obtained from Eqs. (B2), (B3), and 
(B6). The amplitude may be expressed as a function 
of rr and (ry + Tp), and for the larger values of rz it is 
very nearly linear in (ty + rp). Again, the comparison 
with the analog data shows excellent agreement. The 
small inherent delays in the simulation, which give 
rise to some small disagreements in dwell time when 
the theoretical rp = 0, produce no similar discernible 
effect on limit cycle amplitude. 

It may be noted that for finite values of the system 
parameters rz, Tp, and rp, the normalized dwell time 
Ty, and hence the limit-cycle amplitude x, remain finite. 
By letting k — O in the relationships developed it 
appears that the dwell time ¢,, and limit cycle amplitude 
0 become infinite when tz = tp + (tr/2). It can be 
shown that the results of letting k — O are valid for 
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Tr large. 
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k = 0. The paradox that for any k > 0 there is no 
0 there may be a singularity, 
TH /k, so that 


singularity, while for k 
is resolved by the observation that ty 
unless ty approaches zero faster than k, fy will become 
infinite. 


The curve ft; tp + (tr/2) is the only stability 


boundary when k = 0. For values of ¢, greater than 
this, fy > 0; for smaller values of t,, ty < O while near 
the boundary ty > ©. Thus for k = 0, the boundary 
represents a discontinuity rather than a zero. On the 
other hand, when k > 0, the aerodynamic damping 
provides sufficient lead that not only can ¢;, be less than 
tp + (tr/2), but negative ¢, (lag) is allowable and some- 
times desirable. In the analog computer runs with 
k = O, when tz, was less than this required value the 
system always diverged, while for values of ¢, greater 
than the critical value a stable limit cycle always 
occurred. 


The case K < 0 was not investigated analytically but 
was examined on the computer. Such values represent 
aerodynamically unstable vehicles. Even with the 
addition of lead (7, > O) the entire system was un- 
stable, in the sense of diverging to infinity (or, at least, 
increasing until various parts of the simulation satu- 
rated). In some cases, a limit cycle was reached, but 
when sufficiently large step disturbances were intro- 
duced the system always diverged. In every case, the 
step required to cause divergence was less than one x- 
unit. 


With positive aerodynamic damping (K > 0) a limit 
cycle was reached in every case in which the simulation 
was satisfactory. Furthermore, this limit cycle ap- 
peared to be dynamically stable in every case; that is, 
the system recovered from a step. The step disturb- 
ances applied with K > 0 werex = 1. For the illustra- 
tive example discussed below, L = 750 ft.lbs., 7 = 1000 
slug ft.’, so that when k = 1, this is equivalent to about 
a 40° roll step. Furthermore, these steps were usually 
applied at the maximum value of x. This forcefully 
demonstrates the importance of consideration of the 
aerodynamic damping. The data indicate that the 
recovery from a step is dependent on limit cycle ampli- 
tude: when the limit cycle amplitude is relatively 
large, recovery is rapid. 
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Phase-Plane Trajectories 


Although the primary analytic tool for this study was 
a time domain approach, additional insight into the 
contribution of the various system parameters on the 
system behavior is given by a series of illustrative 
phase-plane trajectories. 

In Fig. 9 the effect of lead and delay on the system 
stability is shown for the case of no aerodynamic damp- 
ing. In this case, a trajectory is composed of a sequence 
of segments derived from the parabola 


Le 6° 
—_ + constant 
and its mirror image. When t, = tp = tr = O, one 


obtains a vortex, that is, a limit cycle whose amplitude 
and frequency depend on the initial conditions. For 
lead alone, the origin is a stable limit point, while for 
delay alone an unstable oscillation results. The de- 
stabilizing effect of delay may be compensated by 
sufficient lead, tz > tp + (tr/2), resulting in a stable 
limit cycie of finite amplitude. 

The effect of aerodynamic damping on the system 
stability is similar to that of lead, as is clearly shown in 
Fig. 10. The trajectory for the case of positive aero- 
dynamic damping with rr = 0 is generated from seg- 
ments of the following curve obtained from Eq. (B1): 


gn us’ th i + x’| + constant 


Figure 11 shows the effect of increasing lead. When 
all control system parameters are zero, the origin is a 
stable limit point. With positive lead (rz, > 0) the 
origin remains a limit point; however, the approach 
to this point is now along the line 


x + 7rzx’ = 0 


with the control flip-flopping at an infinite frequency 
with zero dwell time, but nonzero net torque. Actually, 
the control system characteristic rz, uniquely deter- 
mines the rate of approach to the limit point at the 
origin. In fact, once the control reversal line 


x + ryx’ = 0 
is reached, say at x = x», the system behavior is charac- 
terized by 

signa” 


or, in real time and angle, 
—t 
6/0) =e th 


so that the total system response from this instant on is 
independent of the magnitude (but not the existence) 
of the aerodynamic damping. 

When the lead time becomes infinite, the system is 
controlled by velocity only and every point of the line 
x’ = Oisa limit point. The actual limit point reached 
by the system depends only on the initial conditions. 
Upon reaching this limit point, the control flip-flops at 
infinite frequency with zero dwell time and zero net 
torque. 


Figure 12 indicates the effect of system delay on the 
system behavior. With delay, one always obtains a 
limit cycle rather than the limit points shown in Fig. 11. 
Comparison with Fig. 9 shows the stabilizing effect of 
aerodynamic damping; whenever damping is present a 
stable limit cycle will be obtained. Of course, with 















































Fic.9. Effect of lead and system delay; no aerodynamic damping 
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Fic 10. Effect of aerodynamic damping. 


little aerodynamic damping and no lead, the magnitude 
and period of the stable limit cycle might be so large as 
to make the system physically unworkable. Although 
not presented in Fig. 12, the result of system delay 
with an infinite lead time is a limit cycle of approximate 
amplitude 0.25, about an arbitrary point on the x-axis 
whose value depends on the initial conditions. 

No phase-plane plots have been presented for the 
condition tg > 0, since the qualitative effect in the 
phase plane of such values of rz is similar to that of a 
delay, so that the tedium of obtaining quantitative 
results did not appear to be warranted. 


Application to Flight 


One of the more promising vehicle types for hyper- 
sonic flight, the boost-glide vehicle, was examined to 
underscore the importance of considering aerodynamic 
damping in the design of practical flip-flop control 
systems. 

The boost-glide vehicle, as its name implies, is hurtled 
to very high altitude and hypersonic speed during a 
brief initial boost phase. It then glides along a flat 
glide trajectory toward earth, utilizing lift to stay well 
within but near the upper altitude bound of the so-called 
“corridor of flight,’’ thus avoiding excessive aero- 
dynamic heating. As a vehicle type, it represents not 
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TABLE | 
h, ft. k, sec.~! 
225,000 0.024 
181,000 0.12 
150,000 0.4 
117,000 2.0 


101,000 4.0 


only an efficient flying machine in its own right, but its 
trajectory is also representative of what might be ex- 
pected for re-entering satellite and space vehicles, either 
manned or unmanned. The speeds, flight times and 
altitudes encountered by a hypothetical boost-glide 
vehicle are indicated in Fig. 5 of reference 4. 

At the altitudes and speeds involved in this trajec- 
tory, particularly the higher altitude end, it is more 
than likely that reaction controls, rather than con- 
ventional ailerons, will supply the roll control torque 
and roll stabilization. 

The damping-in-roll-derivative for a hypothetical 
boost-glide vehicle, shown in Fig. 13, was estimated 
from data given in reference 5. When combined with 
the trajectory data of reference 4, and the following 
assumed values for the wing area (S = 460 sq. ft.), wing 
span (6 = 40 ft.) and moment of inertia in roll (J = 1000 
slug ft.”), the values shown in Table | for the normalized 
damping parameter k are obtained. 

For reasonable manned flight a desired cycle time of 
the flip-flop roll-control system would be about 2.5 sec. 
A hold time of about 1 sec. would therefore be required, 
if the realistic rise and delay times of 0.1 sec. were 
assumed as characteristic of the system. The designer 
now may follow one of two alternatives: 

(a) He can disregard the damping and design the 
system lead time for k = 0, allowing the control dwell 
time to vary along the trajectory; or 

(b) he can consider the damping and try to hold the 
desired control dwell time by making the lead time vary 
along the trajectory. 

Which alternative is chosen will depend strongly on 
the impact that damping has on the dwell time if the 
lead time is held fixed, as well as on the requirement to 
hold the desired cycle time between narrow limits. 

Although the cycle time of the system may appear at 
first glance to be an arbitrary choice left to the design- 
er’s discretion, considerations of frequency coupling 
with other rotational modes (such as pitch or yaw), 
resonance and vibrational damage to internal compo- 
nents, pilot and passenger comfort, ability to use the air- 
frame as a gun, bombing, or rocket-launching platform, 
and even some forms of buffet and aeroelastic behavior 
prescribe rather narrow bounds on the designer's 
allowable choice of roll frequency and amplitude. 

Utilizing the equations developed in Appendices A 
and B, the system was analyzed for both alternatives. 
The variation of dwell time and rolling amplitude for 
both conditions is summarized graphically in Fig. 14. 
The required variation of lead time to hold the dwell 
time at the desired value of 1 sec. is given in Table 2. 

The effect of damping on the system performance is 
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Fic. 12. Effect of system delay. 
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Fic. 14. Effect of aerodynamic damping on hypothetical system 


performance. 


strong enough that it cannot be neglected. As might 
be expected, its effect is akin to that of an added lead 
time; thus, increasing the aerodynamic damping re- 
quires decreasing the lead, or even adding a lag to the 
system to maintain the desired characteristics. 

It is apparent from Fig. 14 that designing the system 
with a fixed lead time—the simplest procedure—and 
letting the dwell time vary results in a system that can 
very easily give trouble in flight, particularly when 
frequency coupling with other rotational modes, vibra- 
tion, and manned utility of the vehicle are considered. 

Designing the system to hold dwell time constant 
also appears to be a formidable problem. It is apparent 
that some form of variable gain in the lead circuit will be 
required. The variation of gain is so great, that cou- 
pled with the uncertainties of performance, it would 
appear that to utilize flip-flop control systems properly 
will require that adaptive systems rather than pro- 
grammed variable-gain systems be considered. 

It is interesting to note the effect of damping on the 
limit cycle amplitude. As can be seen, given constant 
rise, delay, and dwell times, damping eats into the roll- 
ing torque in a substantial manner, reducing it approxi- 
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TABLE 2 
k, sec.—! t, sec. 
0 0.174 
0.024 0.174 
0.120 0.154 
0.4 0.108 
2.0 —0.150 


4.0 —0.214 
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mately 20% per unit of the normalized parameter . 
Its qualitative effect, if lead time is held constant and 
dwell time is allowed to vary, can be intuitively inferred. 
The magnitude of the effect is rather surprising, how- 
ever, indicating that long recovery times from a 
transient may be expected from a system with constant 
lead time. 

Several other possible combinations of lead, rise, and 
dwell times were also investigated as part of this study. 
All exhibited the behavior given in this representative 
example. It is interesting to note, however, that when 
the normalizing parameter & is less than 0.075, the 
damping need not be considered, and the system be- 
haves as if there were no aerodynamic damping. 


Conclusions and Recommendations 


Closed-form analytical solutions have been found 
which relate the various parameters in closed-loop flip- 
flop control systems. The relationships among lead 
time, rise time, delay time, dwell time, and limit-cycle 
amplitude have been determined as functions of inertia, 
aerodynamic damping, and available control moment. 

A strong dependence of system performance on aero- 
dynamic damping has been demonstrated. Damping 
is equivalent to a lead time, and when heavy aero- 
dynamic damping is present a lag in the feedback loop 
may be necessary for desired system performance. 

For those vehicles operating over wide ranges of speed 
and altitude, optimal control systems appear possible 
only when adaptive systems are used in conjunction 
with the flip-flop. 

Several facets of the problem not investigated in de- 
tail in this study warrant further investigation. Among 
criteria for the dynamic stability of the 

the coexistent limit 
cycles when the small; the 
effects of dead zones and hysteresis on the system per- 


these are: 


limit cycle; stability of two 


control rise time is 
formance; and transient response in recovering from a 


disturbance, including disturbances other than steps. 


Appendix A 


The relationships among the time constants /;, tp, fr, 
and ty will be obtained. It has been shown earlier that 
after normalization one obtains the piecewise linear 


system satisfying at various times one of the following: 


x” +x’ = -] (Al) 

x” +x’ = —1 + (27/rpR) (A2) 
oe = (A3) 

x” +x! = 1 — (27/TpR) (A4) 


The general solution of each of these segments may be 
obtained quite simply. 

Assume that a limit cycle exists.* 
the instant that the control in the hard-over negative 


Consider 7 = Oas 


* In this and the following Appendix, necessary conditions for 
a limit cycle are derived, without any investigation of their 
sufficiency. 
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Then for 0 < 7 < 7p, 
Moreover, at 7 OQ, 


position receives the flip signal. 
the system satisfies Eq. (A1). 


x+Tx =0 


Solving (Al) with this initial condition leaves one un- 
determined coefficient. This will be eliminated subse- 
quently, essentially by identifying the initial conditions 
with the final conditions just prior to receiving the flip 
signal. 


In terms of this as yet undetermined coefficient, 


x 77(B + 1) — B+ Be — T ~ 
: & J é . { (AS) 
x’ = —Be'’ -— 1 O<r<rp f 
and for r+ 0, 
x 7,(B + 1) 
’ . : , t (A6) 
x —(B+1)) 
After r Tp the system is described by (A2). The 


initial conditions (at 7 Tp) are obtained from (A5). 
It is convenient to take as a new time origin the instant 


the control starts to move. Then, for 0 < 7 < Tp, 


») 
\ T11(B + 1) — B — tp + + 
‘R 
ad 2 : 2 T 
Be-T — e7—{1+ 7+ (A7) 
TR TR TR 
2 2 2r 
x" — | Be-T? — iad RE + 
TR TR TR 
After r = Tp the system is described by (A3) with 
initial conditions given by (A7) at r = rr. Shifting 
the time origin to the instant the control reaches the +1 
position, one obtains at 7 = Ty, 
vr = 7,(B+1)-—-B—tr—-—2+ 
9 9 
Be"? — — }e*® + — b-** + sq 
TR TR ) (AS) 
>) 9) 
y= — Be-7? — e77TR + e7TH + | 
TR TR. 
At this instant the flip signal occurs again. The 


system is described by (A3) and the initial conditions 
satisfy 

x + 7r,x’ = 0 
It is simpler not to use (AS) as initial conditions but 
instead to carry another constant, B;, also to be deter- 


mined later. In terms of this constant the initial 
conditions are 
x = (Bi -—1)r 
s ah ut (AQ) 
x = —(B, ail 1)f 


Proceeding as before, the conditions at the next trip 


signal are obtained as 


(A10) 


LIP-FLOP 
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Identifying (A10) with (AG), and (A9) with (AS), one 
obtains four linear equations in the two unknowns B 


and 5;. Collecting terms: 


rr — 1 + e(tot+tR+TH))B — 7B, 
9 
—2T; 7 Tea J 2 — r Tr ] — <¢ TR) — TH 
‘R 
9 
e-(rot+re+1H)B — By = — — e-TH(1 — e-7) 
‘R 
T1rB— (rp — lL +e DTTRTTH)) B, 
» 
—2r, “T Tp T z —_ ( TE l —'<¢ TR) — TH 
‘RR 
» 
B-—. TOT TE TH) B, —_ e~TH(]| — e-TR) 
TR 


Clearly, the first equation is a linear combination of 
the other three and may be discarded. For the result- 
ing three equations in two unknowns to have a solution 
it is necessary that the determinant of coefficients be 
equal to zero. After some manipulation this condition 


becomes: 


2(1 — T) = 
Tete — To)(l + e OT RT Ta) 411 
(2 ) 
tTrR(1 + eSTOtTTRTTH)) — De-TH(] — ¢ R 


An explicit solution for ry is obviously impossible. 
However, numerical or graphical solutions for ry can be 
obtained if required; for example, see Figs. 4 and 5. 

Examination of Eq. (All) shows that ry becomes 
infinite only when 7, > ©. The parameter 7, becomes 
infinite when 


( 4 rer ‘ y > 
Tr(1 + e TDT+TR TH )— 2¢e TH(] —e TR) — () 


For non-negative values of rp, rr, and ry the left side of 
this expression exceeds 


Tr — 211 —e 7) > 0O whenever rp> 1.59 
Since the numerator of the right side of (A1l1) generally 
does not go through zero when the denominator does, 
it follows that 7, goes from + © to — © at this point. 
This is illustrated in Fig. 5 of the text. 

The “‘stability’’ boundary of the system is obtained 
by setting ry = Oin (A11), obtaining 
TRTTD ') 


—TrTp(l +e 


Tr(l + eTRTTD)) — 2(] — e-TR) 


2(1— rz) = (A12) 


which is shown in Fig. 6 of the text, and which exhibits 
properties similar to (A11). 


Special Cases 
(a) tre = 0. The derivation of (All) remains valid 
when 7z = 0, so that one obtains by passing to the limit 


in (All): 


es (TH Tp)(1 +e \TDTTH)) 
2(1 = TL) = > lien dveurt 9 i ’ fs, = 0 
l +e TD H — Ve H 
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Fic. 15. 


When k = 0 the form of the Eqs. (A1l)- 
Proceeding in a manner similar to that 


(6) k = 0. 
(A4) changes. 
used in obtaining (A1l1) we obtain 
= tr + dt rtp aE Glut k ae 


3(2t, — 2tp — tr) 


ty (A13) 
It is interesting to note that by substituting + = At in 
(A11) and letting k — 0, one also obtains (A13). The 
qualitative differences between (All) and (A13) have 
been discussed in the text. 


Determination of B and B, 

For use in Appendix B, it will be convenient to have 
solutions for B and B,. Of the four equations connect- 
ing B and B, presented above, the second and fourth do 
not contain r;. These two equations are repeated here. 


9 


Be-(tot+tr+TH) — By = — — e-TH(] — eR) 
TR 
9 
B — Bye-tot+tR+1TH) = — — e-TH(] — eR) 
TR 
Clearly: 
2 i-—<¢" , 
—B = B, = —e-™- e*) ___ (414) 
TR 1 + e-(tottRt+TH) 
The fact that B; = —B verifies the physically obvious 


fact that the limit cycle is symmetric about x = 0; that 
is, there is no bias in the trim position. 


Appendix B 


The equations defining the limit cycle amplitude will 
be derived. The equivalent of the phase plane (x,x’) 
trajectory of the solution is used. The phase-plane 
trajectory might be obtained in the usual fashion by 
making the change of variables 


f dx 
i = 
2X 4 


in the Eqs. (Al)-(A4). No difficulties are en- 


countered for the segments (Al) and (A3), but in the 
segments (A2) and (A4) a nonlinear differential equa- 
tion is obtained; however, this equation can be solved. 
A less erudite but in this case simpler method is to 
eliminate 7 from the explicit solutions (A5) and (A7) 
after substituting for B from (A14). 


NOVEMBER, 1960 


SCIENCES- 


Then, as long as the normalized control torque is — 1, 


' Tot+T 
x = —x’ + In|1 + x'| + “~~ + 
tnl+e (ro+TR+TH) 
in| = | (B1) 
2 L=—¢ = 


Thus, if a maximum or minimum occurs on this segment 
of the control cycle, x’ = 0 and 


l+e (TO+TRt+TH) 
es Tp + Te + - E | (B2) 


2 2 l1—e TR 


The time variable 7 does not appear in the phase plane. 
Hence, one cannot determine in the phase plane if a 
maximum or minimum does occur while (B1) is appli- 
cable. From (A5), however, it is clear that at most one 
stationary point (at r = +In(—8B)) can occur while 
(A5) is valid. Since (A5) is valid as long as the control 
torque is —1(—rTy < 7 < rp), there can be at most one 
stationary point on this segment occurring when 


=e < 7 = In (-—B) = 


2 1—e TR 
In —T 
Tr 1+ ettotrRt+rTH) ad 


During the control torque rise from —1 to +1, Eqs. 
(A7) are valid. Adding these two equations and sub- 
stituting for B, one obtains a quadratic equation for r. 
Solving, one finds that when x’ = 0, 


9 
TR TR TRTH TRTD 
f=— + + + —— + rt 
4 2 2 


TD (B3) 


lA 


= St y (B4) 
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Substituting in the second equation of (A7) with x’ = 0, 


»—(tH+T 
1 = yje™” = Les : “0 e-Tk/2 = q _ (B5) 
1 + e-(ta+rbp+TR) 

Clearly, a> 0. It may easily be shown that a < 1. 
The functional relation between a and y is shown in 
Fig. 15. Actually there is a second solution y < 0 to 
this equation. This second solution may be shown to 
satisfv ly| > tr/2, so that from (B4) it corresponds to 
an inadmissable value of r. 

The positive value of y may be greater than 7r,/2. 
In this case, no stationary point occurs on the ramp. 
When 0 < y < 7rp/2 it would appear from (B4) that 
two stationary points might occur at both (rz/2) + y 
and (rr/2) — y. However, from Eq. (B4), 

2 
ds ace ‘rR  *H + Tp - 7 (B6) 
+ 2 TR 
so that the value of x is the same at both stationary 
points on the same ramp. This is impossible; for if 
both points were maxima (minima) there would have 
to be a minimum (maximum) in between, while if one 
point is a maximum and the other a minimum their 
values could not be equal. Thus, there can be at most 
one stationary point on the ramp. 

Let us review our status. We have shown that there 
is at most one stationary point on the ramp and at most 
one stationary point on the flat. With B = —A,, the 
explicit solutions of Appendix A show that if a maxi- 
mum occurs on one flat segment a minimum must 
occur on the next flat segment. Between a maximum 
and minimum there must be an even number of maxima 
plus minima. Therefore, there cannot be one station- 
ary point on the ramp. Thus, we have shown that 


there is exactly one stationary point per half cycle. 
If (B3) is satisfied it occurs on the flat segment; if (B3) 
is violated it occurs on the ramp. It is necessary to 
determine whether it occurs on the flat or ramp segment, 
since different formulae for the limit cycle amplitude 
must be applied in each case. 

Figure 16 is a graphic portrayal of (B3), the condition 
determining where the maximum occurs. Thus to find 
the limit cycle amplitude one first refers to Fig. 16. 
If the peaks occur on the flat segments, the limit cycle 
amplitude is given by Eq. (B2). If the peaks occur on 
the ramps, one evaluates a and uses either Fig. 15 or 
Eq. (B5) to obtain y._ [For normal purposes, the graph 
is sufficiently accurate so that numerical solution of 
(B5) is unnecessary.| After obtaining y, substitution 
in Eq. (B6) gives the limit cycle amplitude. Figures 7 
and 8 presented in the text include points falling in 
both regions—peaks on the flats and peaks on the 
ramps—and no discontinuity occurs at the junction. 
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Elastic Instability of a Cylindrical Shell 
Under Arbitrary Circumferential 
Variation of Axial Stress 


P. P. BISLAARD* ano R. H. GALLAGHER** 
Bell Aderosystems Company 


Summary 


The buckling stress of a cylindrical shell is determined for 
various circumferential distributions of axial stress. The de- 
velopment is based on small deflection theory, with finite dif- 
ference techniques being applied to derive the general form of 
the governing algebraic linear homogeneous equations. The 
equations are cast in both the determinant and matrix forms, 
the latter being suitable for the commonly used matrix iteration 
solution. Numerical results are obtained by means of a high- 
speed digital computer for cases wherein the stress distribution 
is described by 32 circumferential elements. It is found that 
the buckling stress for pure bending, and for more complicated 
circumferential stress distributions of the axial stress as well, 
is not much higher than for uniform axial compression. Cases 
are indicated where small deflection theory is directly applicable. 


Introduction 


A“ INTENSIFIED INTEREST in the behavior under 
load of structures composed of cylindrical shell 
elements has been brought about by the configurations 
featured by missiles and hypersonic manned aircraft. 
The circular cross-sectioned fuselage is a predominant 
part of both types of vehicles, and to make possible 
satisfactory performance it is necessary to proportion 
the structural elements of the fuselage in an efficient 
manner. For elements about which there is little 
known concerning strength and deformation character- 
istics, an efficient design can be evolved only when reli- 
able information is made available. 

The unreinforced, thin-walled cylindrical shell is of 
particular concern in this regard. At present, to de- 
sign this type of structure for the comparatively simple 
case of uniform axial compression one must rely mainly 
on test results, since no suitable theory is available 
for the direct computation of the critical stress. Fur- 
thermore, the test results extend only to the case of 
pure bending, with scant data pertaining to intermedi- 
ate conditions of combined bending and axial com- 
pression. Yet, even if this information were available 
it would not be sufficient for conditions of axial stress 
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brought about by the circumferential temperature 
gradients due to aerodynamic heating. Such gra- 
dients can arise as a result of differences between the 
upper and lower surface heating rates as caused by a 
positive angle of attack, or as a result of the heating 
conditions that prevail in a partly-full integral fuselage 
propellant tank. In either case, there are developed 
significant axial thermal stresses that are nonuni- 
formly distributed around the circumference. In view 
of the arbitrary nature of thermal stress profiles, as 
well as the changes that take place in the total stress 
profile as external load stresses are superimposed, 
linear stress distribution data would not suffice for 
the determination of critical stresses which may be 
supposed to depend on the nature of the stress dis- 
tribution. 

Rather than to try to solve the problem of deter- 
mining the critical stress for a circumferentially vary- 
ing axial stress system by means of large deflection 
theory, as a first approach small deflection theory is 
used. This will give some indication of the significance 
of arbitrarily distributed axial stress as compared to 
homogeneous axial stress. Also, there are at least two 
important situations where the usefulness of the small 
deflection approach is directly apparent. The first 
situation relates to pressurized cylindrical shells. It 
has been contended, and to a certain extent verified,! 
that when a given amount of internal pressure is pres- 
ent, a cylindrical shell under uniform axial compression 
will buckle at the stress level predicted by small de- 
flection theory. Greater pressurization will not yield 
an improved critical stress. This is in accordance 
with small deflection theory which does not indicate 
any greater critical stress due to internal pressure. 

A second viewpoint associated with the usefulness 
of the general approach used in this paper is derivable 
from the work in progress on the plastic buckling of 
cylindrical shells. Here, due to the relative ‘‘thick- 
ness” of some shells, it is found that small deflection 
theory, combined with the theory of plastic buckling of 
shells, leads to results in accordance with experi- 
ments.”~* Also, for thinner shells, small deflection 
theory leads to a plastic reduction factor that can be 
applied to the elastic buckling stress from large de- 
flection theory.® 

The method of solution consists of first a reduction 
of Donnell’s partial differential equation for the radial 
deflection w(x, y) to an ordinary differential equation 
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for the circumferentially-dependent radial deflection 
S(y), and then the application of second-order finite 
differences to obtain general expressions in an algebraic 
form. The latter technique is facilitated by the fact 
that the shell is ‘‘round,”’ whereby the need to specify 
“outer” points or to establish boundary conditions is 
obviated. By dividing the into 2n 
spacings the approach leads to an (m + 1)th order 
determinant for stress distributions symmetrical to the 
Z-axis (Fig. 1). This determinant, when equated to 
zero, yields the buckling stress. However, since high- 
speed digital computers are more amenable to solution 
procedures that employ matrix iteration, the matrix 
form of the governing equations is indicated. 

In accomplishing the solution to numerical prob- 
lems, the authors have chosen to embody the final re- 
sult in the coefficient AK which is defined in conjunction 
with its use in the equation 


circumference 


Gn = Kitt/a) (1) 


This is the most common form of the thin-walled 
cylinder design equation, and when used in conjunc- 
tion with test data it is found that A has a dependence 
on the thickness-to-radius ratio. The dependence is 
not significantly represented in results that are based 


on small deflection theory. 


Derivations 


Assuming a loading by axial stresses o, which are 
distributed arbitrarily in the circumferential direction 
but are assumed to be constant in the axial direction 


(see Fig. 1), Donnell’s equation becomes 


i | = — 2] Ou! + ft wr 07x — 
Ww o = (Z 
at? Ox dD * Or? 
where 

Vi = [(07/dx?) + (0?/0s”)]? 

vs _ Viv 

D = shell flexural rigidity = Et*/12(t — v’) 

a = shell radius 

t = wall thickness 


Eq. (2) can be reduced to an ordinary differential 
equation if an assumption is first made concerning the 


longitudinal shape of the radial deflection w. For 
this purpose it is assumed that 
w= Ssin Ax (3) 


where 


A = az/l 

1 = the half wavelength of the buckles in the longi- 
tudinal direction 

S = the circumferential buckle shape function (a 


function onty of the circumferential coordi- 
nate, s) 


) 


Substitution of Eq. (3) into Eq. (2) results in 


8S ane a°s + 6a! d's 
; " ds? ' ds4 
aS fS _12(1 — »*) . 
4)" : age SAt + 
ds® ds* a’*t? 
1*( pS) 14( pS) 
° | ~asn +ai—— _ yi eo @ 
ds? ds 
where p 2 O./Ouec (5) 
@ = ig,,,,/D (6) 


p is therefore a function that describes the circum- 
ferential variation of axial stress and is obtained by 
normalizing all axial stresses on the maximum axial 
Stress, Omar. Since p and S are functions of s only, the 
product pS is also a function of s only and will be de- 


scribed as Z, whence 


Z = f(s) = pS (7) 


Utilizing the superscripts II, IV, VI, and VIII to indi- 
cate the second, fourth, sixth, and eighth derivatives 
with respect to s, Eq. (4) then becomes 


sv™ — ars™ + acs’? — aes" + 


12(1 — v?)A4 — 7 
s+ — — #°\%p, | S — 


at 
®°Z'Y + 26°17" = 0 (8) 
In order to write this equation in terms of finite dif- 


ferences, the circumference is divided into 2” equal 
spacings h (= ra/n), as illustrated in Fig. 2. To make 
the calculation as accurate as possible, as in an earlier 
The 
derivation of the finite difference expressions for cer- 
tain derivatives can be found in a number of texts.* ® 
Since the expressions for the second and fourth deriva- 


paper’ central second-order differences are used. 


tives are given in reference 9, only those for the sixth 
and eighth derivatives, which occur for S, will be given 


here. They are 


4n8SY' = [—Sy-14 + 12 Se; — 52 Se-2 + 
116 S.—, — 150 S: + 116 Sus: = 
52 Sete + 12 Sere — Seta] 

[—S-s + 13 Ss — 69 S-3 + 
204 S,-2 — 378 S,-1 + 462 X 
S, — 378 Seti + 204 Stas — 
69 Sirs + 13 Seas — Servs] 


9),8 eVIII 
3h8S (9) 


With the expression for the second and fourth deriva- 
tives as given in reference 9 and taking account of Eq. 
(7), from which 
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Zn = PnSn (10) 


the second and fourth derivatives of the variable Z are 
determined as 


1272" ad [ — pr—-2Sk-2 + 16 px—1Sx—-1 = 
B0 puSe + 16 pet iSet1 — peteSe+2] 
GHZ = [—py—2Se—3 + 12pp-2S-2 — > (11) 


39 px- We-1 + 56 pp Sy —_ 39 pet Set 
+ 12pptoSet2 — pet aSets | 
Insertion of all these expressions, including those for 
S" and S'Y, into Eq. (8) yields 


—4[Sps + Sersl + CilSea + Seta] — 
CoS-3 + CoSh-2 — CrSi-1 + 


CsS; —- CeSe4s + CrSit+2 — CaSx4 3 = 0 (12) 
where 
CG; = 52 + 122 
C2 = 276 + 144y + 12? —_ 211 p.x-3 
C; = 816 + 624u + 144u? + 4y3 — 
(u + 12) 2Mp,.—2 
C, = 1512 + 13924 + 468u? + 64yn? — 
(16u + 39) 21 p..-1 
C; = 1848 + 1800u + 672y? + 120u* + (13) 
12u4 + 12u0 — (3u2 + 15 + 28)4IIp,, 
Cs = 1512 + 1392u + 468yu? + 64? — 
(162 + 39)21p,x41 
C; = 816 + 624u + 144? + 4p? — 
(u + 12)2TI p,x+2 
Cs = 276 + 1444 + 12u? — 2 psp42 
and 
p = Nh (14) 
08 = 12[r4(1 — v?)/n4](a/t)? (15) 
II = A2h*h? (16) 


An applicable form of Eq. (12) can be written for each 
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finite difference element, resulting in a set of 2m linear 
homogeneous equations in the unknown S's. In order 
to find nonzero solutions for the function S the deter- 
minant of their coefficients must be zero. Hence, 
expansion of the determinant yields a polynomial ex- 
pression with II (the critical stress parameter) as the 
unknown, and solutions for II can be attempted after 
the condition that the expression must equal zero has 
been applied. 

If attention is directed circumferential 
stress variations that are symmetrical about the Z-axis 
(Fig. 1), the most commonly encountered condition, 
the general form of the determinant array can easily 
be written. Due to symmetry, only the z + 1 elements 
of the half circle need be considered and the resulting 
form of the array is as shown in Fig. 3. 

With the development of high-speed digital compu- 
ters, it is feasible to effect the numerical solution to a 
problem of the type represented by Eqs. (12) for large 
values of nm. The most efficient solution procedure for 
digital computers is often that which employs matrix 
algebra and matrix iteration. The transformation of 


toward 
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ELASTIC INSTABILITY 
Eq. (12) into a statement in matrix form is accom- 
plished as follows. Denoting an individual multiplier 
of S in Eq. (12) as &;;, it can be seen from Eq. (13) 
that any such term can be expressed as: 


Ey = Ay — B,ll (17) 
By use of this form for each term of (12) it is possible 
to write the matrix statement of the total set of Eqs. 


(12) in the form 


({A] — II[B])}S} = 0 (18) 


where [A] and [B] are matrices populated with the 
terms A,; and B,,;, respectively, and {s} is a column 
matrix of the unknown amplitudes S;,..., ¢ we «OA 
description of the iterative procedures that can be used 
for the solution of (18) can be found in reference 12. 

In review of the data necessary for a numerical solu- 
tion, it is seen from (13), (14), (15), and (17) that the 
A,; and B,; terms are derived from geometrical data 
(given a/t, assumed //a) and the stress distribution. 
These terms are therefore constants for a given problem. 
When a solution for II has been accomplished, the criti- 
cal stress Omar follows from Eqs. (6) and (16), so that 
with h = za/n 

Et?/? 


mar ~~ oe II 
” 12(1 — v?)e2h! 


E n* ft \* fi \* 
ie II (19) 
12(1 — v?) r® \a a 


Substitution of the solution for II as well as the value 
of //a and t/a (the input upon which the solution had 
been based) will provide for the expression of omg; in 
terms of the modulus of elasticity. From Eqs. (1) and 
(19), K follows as 


4 t l 2 
- 7 1 = . (EN) . : 
12(1 — v?)r® \a/\a 
. rw 
9.533 XK 10-5 n*([- }| -] IE (20) 
a/\a 


Numerical Results 


Critical stresses were determined by means of the 
above approach for a variety of circumferentially 
varying axial stress systems. Since it is at first neces- 
sary to assume a particular value for the radius to thick- 
ness ratio, a value of a/t = 1000/+/12 was used in most 
cases. This ratio is identical to the one used by Fliigge 
in reference 10. A limited number of determinations 
were made for a/t ratios of 100 and 1000 for the purpose 
of studying the influence of this geometric parameter 
on the values obtained for K. 

To facilitate the numerical computations a program 
was written for the IBM 704 digital computer. The 
solution approach was one of direct matrix iteration. 
The governing matrices of each problem were formu- 
lated by the computer on the basis of the following in- 
put data: 


OF A CYLINDRICAL 


SHELL 857 


(a) An assumption of the buckle length-to-radius 
ratio //a; 

(b) A single geometric parameter 0 
the a/t ratio of the cylinder, 

(c) The m + 1 values of p, 
of the half-section. This describes completely, in non- 
dimensional form, the circumferential stress distribu- 


a function of 


a value for each element 


tion. 

A number of //a ratios must be investigated to ob- 
tain the value of critical stress for any given profile, 
since, in accordance with theoretical assumptions, the 
cylinder is “‘long’’ and is able to buckle at whatever 
ratio is consistent with the smallest critical stress. 
In view of this factor, and in order to provide the maxi- 
mum of useful results and conclusions, the studies were 
directed toward “harmonic” circumferentially varying 


axial stresses. These are defined as 


o; = [cos (ms/a)] Omex (21) 
so that, from Eq. (5), 
p = cos (ms/a) (22) 


The assumed division of the circumference into 32 
elements will make the mesh closely representative of 
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the cos ¢ = cos (ms/a) (m = 1) harmonic, which de- 


fines a stress distribution due to pure flexure. Pro- 
files that for data presentation purposes are given by 
cos 8p = cos (8s/a) (m = 8) or cos 16¢ = cos (16s/a) 
(m = 16) are truly representative only of the distri- 
butions shown in Fig. 4. 

As a check case, the critical stresses of a cylinder 
under uniform axial compression were determined. 
Values obtained for nine assumed //a ratios have been 
plotted in Fig. 5. For comparison purposes, the solu- 
tion based on Eq. (264) of reference 11 is also shown. 
Taking account of the results shown in Fig. 240 of 
reference 11, the present solution can be considered as 
satisfactory. 

The investigation of the cos (s/a) = cos @ (bending) 
case provided results of considerable interest, since 
certain opinions exist relative to this profile despite 
the absence of a comprehensive solution. For instance, 
Fliigge™ studied this case for an assumed //a ratio of 
mw and obtained a critical stress at the crown, Omaz; 
which was 1.3 times the cylinder critical stress under 
uniform axial compression. 

Seven assumed axial buckle lengths were studied for 
the case of bending loads, and the results have been 
plotted in Fig. 5. It is seen that a correspondence 
exists between the result of Fliigge and of the present 
procedure. For lower values of assumed //a ratios, 
however, only a small margin exists between the bend- 
ing and uniform compression critical stresses; for an 
assumed //a of 0.8, the critical stress due to bending is 
found to be only 4% greater than the smallest uniform 
axial compression critical stress. Hence, within the 
framework of linear theory, buckling under either flex- 
ure or direct stress occurs at essentially the same ex- 
treme fiber stress level. It is of interest to note that 


in a recently published paper by Abir and Nardo™ 
essentially the same conclusion was drawn based on a 
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series solution to Donnell’s equation. Abir and Nardo 
obtained solutions for linear variations in stress such 
as would result from the simultaneous application of a 
longitudinal force and bending moment to the cylinder. 
The applicability of this conclusion to realistic situa- 
tions was discussed in our Introduction. 

Solutions were also obtained for cos (2s/a) = cos 2¢ 
(m = 2), the assumed values of //a being the same as 
those chosen for the case of bending. These resuits 
are also shown in Fig. 5, and again, it is indicated that 
the critical stress is not greatlv in excess of the value for 
uniform axial compression; the ratio between the mini- 
mum buckling stress and that for uniform compression 
is 1.08. For cos 8¢ and cos 16¢ (see Fig. 4) with only 
two buckling lengths investigated, the minimum criti- 
cal stress ratios were 1.180 and 1.263, respectively. 

Two additional illustrations that are of interest can 
be constructed from the results given in Fig. 5. Fig. 6 
shows the minimum values of K = o,.a/Et plotted 
versus the orders of harmonics, m. These minimum 
values are available only for the following harmonics: 
0, 1, 2, 8, and 16. 

Fig. 7, using the same abscissa (m), is a plot of the 
ratio of the critical stress of the harmonics to that which 
is applicable to uniform axial compression. 

Finally, determinations were made of the critical 
stresses of uniform axial compression for a/t ratios of 
100 and 1000, with an //a ratio of 1.0 being assumed in 
each case. Fig. 8 indicates the results, together with 
the point associated with an a/t ratio of 1000/ y¥ 12 
(= 288). It is indicated that the radius-to-thickness 
ratio has little effect on the critical stress connected 
with this //a ratio. 

Special note must be taken of the results shown herein 
for the case of pure bending. The contention is not 
made that the actual critical stress of long, unpres- 
surized cylinders in bending is essentially the same as 
that for uniform compression; tests have verified that 
the bending strength (with respect to stability) is usually 
from 30 to 40% greater than the uniform compression 
strength. The actual behavior can be _ predicted 
analytically, however, only by means of a nonlinear, 
finite-deflection theory which must assume initial 
imperfections which are not known beforehand. 

It has been established, however, that for conditions 


under which small deflection theory applies, no signi- 
(Continued on page 866) 
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A Simple Method of Matric Structural 


Analysis: Part V—Structures Containing Plate 
Klements of Arbitrary Shape and Thickness 


BERTRAM 


KLEIN* 


Convair, A Division of General Dynamics Corporation 


Summary 
Simple equations are used to treat plate elements. The equa- 
tions are basic equilibrium and force-displacement equations. 
The plate elements may be of any shape and have any thickness 
variation. The stresses in the elements may vary in a complex 
pattern. There need not be any edge members attached to the 
plates. All equations are simple to derive and simple in prin- 
ciple. Most equations contain definite integrals. Since in gen- 
eral the shapes and thicknesses of the plate elements are arbi- 
trary, a simple numerical integration scheme is used to evaluate 
these integrals. Numerous numerical examples are worked out 
to illustrate the mechanics of the method and the accuracy attain- 


able. 
Symbols 

A = stringer cross-sectional area or area 
a,b = dimensions 

E = Young’s modulus 

f = function 

G = shear modulus 

h = beam depth or dimension 

z, = length 

M = bending moment 

n = normal loading 

p = coordinate 

q = shear flow 

R = shear force 

§ = coordinate 

t = plate gage 

u,v, %' = displacements in XY, Y, and Z directions 
x, ¥ = rectangular coordinates 

Y, Y, Z = rectangular coordinate axes 
a = coefficient of thermal expansion 
8 = angular rotation 

1 = shear strain 

6 = deflection normal to stringer 
Ay = temperature change 

0, o = angles 

v = Poisson’s ratio 

a = normal stress 

T = shear stress 

Subscripts 

1,2,... = indices or points 

Cz = clockwise 

( = cover 

i,j,k = indices 

n, S$ = in directions m or s 

p = panel 

str = stringer 
w = web 


Received May 21, 1959. 
* Design Specialist. Formerly, Chief of Structures, Solar Air- 
craft Company. 





x,y, xy = in directions with respect to rectangular axes 
@ = in direction defined by @ 
Superscripts 


= average value 


= modified 


Introduction 


Mc” STRUCTURES contain elements that may be 
designated as plate elements. Such elements 
are allowed to take axial as well as shear stresses, in 
contrast to shear panels that are assumed to take shear 
only.'~* Furthermore, in general these plate elements 
may be of any planform and may have arbitrary thick- 
ness variation. The state of stress and strain in an 
element in general is complex in nature. There need 
not be any edge members attached to the plate ele- 
ments. However, in many cases the plate element is 
surrounded by edge-reinforcing members that serve to 
introduce the loads into the plate element, and therefore 
the presence of such stringers is allowed for in the pres- 
ent analysis. The strain in a stringer is assumed to be 
the same as the strain in the plate element along their 
mutual line of attachment. 

Just as for shear panels considered previously, only 
one force-displacement equation for shear need be 
written for a plate element when adjacent reinforcing 
stringers surround this element. When greater detail 
is desired, the plate element may be subdivided into 
Then there are no stringers along 
Conse- 


subplate elements. 
certain edges of these subplate elements. 
quently, axial force-displacement equations must be 
written for these subplate elements in the directions 
of their unsupported edges. 

Each plate (or subplate) element must be in equilib- 
Three overall equations of equilibrium may be 
Since in general the 


rium. 
written for the general element. 
corner angles of plate elements differ from 90°, i.e., 
from right angles, ‘‘Mohr’s-circle’’ equations of equilib- 
rium are required for infinitesimal elements to relate 
stresses in different directions in the corners of plate 
(or subplate) elements. 

On occasion one may desire to consider the bending 
resistance of the edge stringers to normal loadings act- 
ing on the edges of plate (or subplate) elements. In 
that case one adds additional equations of equilibrium 
and force displacement describing the bending behavior 
of the stringers now acting as beams. 
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Fic. 1. General element geometry and loading. 











The correct inclusion of the effects of thermal load- 
ings on plate (or subplate) elements is simple and di- 
rect. One adds certain integrals to the axial force- 
displacement equations involving the strains along the 
edges of the elements. 

It is convenient to employ a fixed set of rectangular 
axes, X and Y, for locating points on the structure. 
Then certain of the integrals and other quantities 
appearing in the above-mentioned equations may be 
taken as referred to these axes. 

It is important to realize that the method presented 
here is not a finite-difference method. In finite- 
difference methods, one replaces differential operators 
appearing in differential equations by the appropriate 
difference equivalents, and proceeds to satisfy the re- 
sulting equations at finite points. In contrast, in the 
method advocated here the structure is broken down 
into a consistent set of large, interlocking elements. 
Then certain integrals involved in defining the behavior 
of these elements are evaluated. In general, these 
integrals must be evaluated by a numerical procedure 
such as Simpson’s rule using the values of functions 
at discrete points. 


Development 


(D Plate-Element Equations 


(1) Equilibrium—Consider the general plate element 
shown in Fig. 1. The positive senses of the edge load- 
ings are the directions given by the arrows, i.e., tension 
is positive for the normal loading and shear is positive 
when the loading causes clockwise rotation about the 
c.g. of the element. From elementary considera- 
tions we derive the three overall equations of equilib- 
rium: 
vertical (positive in Y direction), 


f cxtondx + fext T dy = 0 (1) 
horizontal (positive in X direction), 
— f oxtondy + fextt dx = 0 (2) 


moment (positive from X to Y), 


f cuton(x dx + ydy) + fextr(x dy — ydx) = 0 (8) 


All the above integrals are line integrals evaluated in a 
clockwise direction. 

In addition to the above equations for overall equilib- 
rium of the plate element, one must consider the equilib- 
rium of differential triangular elements taken at the 
corners of the plate element. The following Mohr’s- 
circle equilibrium equations are written for elements 
such as those shown in Fig. 2: 


Cn\ _ sin? @ cos? @ — sin 2¢ O1 
T? —'!/, sin 2 '/2sin 26 cos 2¢/ | a2 (4) 


Notice that the angle ¢ is defined in two ways in the 
figure. 

(2) Force Displacement 
tion is derived (see appendix) : 


1 
ff - Try dX dy = Ff cr(u dx — v dy) (5) 
Ap G 7 


The integral on the left is the integral over the area of 
the plate element. 7x, is positive as shown in Fig. 1. 
The integral involving displacements is a line integral 
taken positive in the clockwise direction around the 
edges of the element, i.e., in the positive direction of 
the shears along the edges, and u and v are displace- 
ments in the directions of the coordinate axes XY and 


for shear the following equa- 


Y, respectively. 
For axial strains along the edges of a plate element, 
the force-displacement equation is written: 


























Fic. 2. Differential triangles relating stresses. 
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Fic. 3. Edge member geometry and loading. 
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(u; — U;) cos 6 + (v; — v;) sin 8 = 


F] ij 
f e, ds -f [(o; — von)/E]ds (6) 


Here @ is the angle measured counterclockwise from the 
X axis to the vector in Fig. 1 in the direction of the 
edge 1-j, where j is at the head end and 7 is at the tail 
end of the vector. Notice that the above equation 
is common to the plate element and to the supporting 
stringer if one is present. 

If an edge of a plate element is completely fixed, the 
strain in the direction of the edge must vanish. Ata 
point on such an edge, one writes: 


0; — vo, = 0 (6a) 


(II) Edge-Member (Stringer) Equations 
(1) Equilibrium—For equilibrium along the stringer 


(see Fig. 3) the equation is 


i 
—(oneA)i + (ont), + f tr, ds = Q (7) 
1 
Perpendicular to the stringer we have 
j 


R,+R;- g,tds = © 


et 


The moments about point / are 
z 
R;Li, + M; -— M; - f $63 ds -= 0 (9) 


(2) Force Displacement—For axial strain the re- 
quired equation is the same as Eq. (6) because the 
strain in the edge member is assumed to be the same 
as the strain in the plate element. For rotation, the 
equation is 


IM 
B; —- Bi — f El ds = 0 (10) 


and for shearing action, 


j *M 
a ee | ~ ds | [= (11 
j i f f EI ds + B; | ds 0 ) 


The last two equations are derived by integrating the 
Bernoulli equation for the bending of a beam. 

All the above equations are nearly sufficient for set- 
ting up problems involving plate elements. Sometimes 
special additional equations are required. In general, 
two equations of equilibrium are required for each joint 
at the intersection of edge members. The strain at a 
point in a stringer may be set equal to the strain in a 
plate element at that point. A special equation may 
be needed to require continuity of a function at a point 
(see appendix). 


(III) Numerical Evaluation of Integrals 


The integrals appearing in all the above equations 
can be evaluated in general by assuming a parabolic 
form for the integrand function and integrating this 
approximating function. This procedure leads to 
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Fic. 4 (top). Plate without stringers. 
Fic. 4(a) (bottom). Ea AT values for plate subjected to heat. 





Simpson's rule for numerical integration if the parabola 
is of second degree. In special cases, higher-order 
parabolas may be needed. 

Assume three equally spaced points 7,7,k. Let the 
spacing between points be called ZL. Thus L = 
Li; = Lj~. Let f(p) be the function to be integrated. 
Then one may derive the following equations by passing 
a second-degree parabola through the three points and 
integrating the resulting function: 


pj 
f f(p) dp a A (Sf; + Sf; = x) = Lf, i(P) (12) 


Dk 
f f(p) dp = 4 (—fi + Sf; + Sf) = Lhia(p) (13) 
D — 


? 


Dk 
f f(p) dp = : (fi + 4f; + fe) = 2LFi-4(p) (14) 


pj p L? : 
[O [1 a0 ap = Fh + 6, - fo = 
pi 0 Z 


di 
pf(p) dp (15) 


p 


pk ep L? 
f f f(p) dp dp = — (—fr + 10f; + 3h) = 
pjiJ/O 24 


f pf(p) dp (16) 
pk 


The last two expressions in Eqs. (15) and (16) represent 
moment terms taken about points j and k respectively. 
Notice that if point 7 is situated on an axis of symmetry 
or antisymmetry, f; = f,. If point & is on an axis of 
symmetry or antisymmetry, Eq. (12) may be rewritten 
(see appendix). 

The double integral appearing in Eq. (4) can be 
evaluated by applying the appropriate equations of 
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First write the relation 


ff Tzy 0X dy = TA, 
Ap 


(12)—-(14) twice. 


Then 7 = (1/12) (57; + 87; — 7) where 7;, 7;, and 7, 


are average values of 7 taken in the directions defined 
by lines passing through the points 7, /, k. 


Numerical Examples 


(1D Plate Element With No Stringers 

A plate element without any edge-reinforcing mem- 
bers is shown in Fig. 4. Advantage is taken of the 
symmetry of the structure and loading which is a uni- 
Only four subplate elements are used 
Twenty- 


form tension. 
to represent the structure in the analysis. 
eight equations are set up as follows: 
(1) Equilibrium: three each for the four elements 
(2) Force displacement: each for the 
four elements; axial edge 
5-6, 7-9, 8-9, 1-7, 2-8, 4-7 ,5-8. special 


shear—one 
ten for edges 1-3, 2-3, 4-6, 
zero strain in X 
direction, one each at points 6 and 9 at wall. 
The corresponding 28 unknown quantities are: 
(1) 12 axial stresses—six each in Y and X directions 
(2) six average shear stresses 
(3) 10 deflections—six in Y direction and four in X 
direction 
Space will not be taken here to represent all the equa- 
Therefore only typical equations are shown 
Consider, for example, equations involving 
The equations are: 


tions. 

below. 

subplate element 4-1-2-5. 
(1) Equilibrium in Y direction: 

ys) + h, T4_5 acs h, Sapp. 


(hy 3)(Gy2 + 3045 — 


in X direction: 


I 


(hy /12)(5or4 ++ Sor = Ox) + h, To 5 0 


moments about point 5: 


(h,?/120)(3loy + 460,5—17 0) — 
(hy?/24) (Tora + Gor — o76) = (1/2)hz*Oapp. 


(2) Force displacement 


(1 '12) (872 eee 73 6)hzhy 'G > (3h, 5) x 
[(u2 + 2u5) — (uw, + 2u4)] + 
(hy/12)[5(v%, — v5) + S(v2 — v5) ] 


shear: 


edge 6-4: 


9) 


“4 = 


(hy/3E) |(oys — vors) + 4(¢y5 — vor) + 


(ays ad Vore) | 
edge 1-7: 


—U = (h, BE) |—voy + 4( org = VOys) + (O27 — voyz) | 


Y1= (h,/3E) (on + 40,2 + Coys) 


edge 2-8: 


--u2 = (h,/3E)[—voy. + 4(o25 — voy) + (e273 — voys 


1960 


NOVEMBER, 


SCIENCES- 


The solution of the 28 equations is: 


o,2 = 0.27390871 O14 = —0.0016921038 
0,3 = 0.34517339 O25 = 0.0080519664 
o,, = 0.28555564 o:x6 = 0.093772931 
o,6 = 0.28131880 O27 = —0.021083315 


0.022147683 
0.081279034 


os = 0.29815445 a3 = 
Cy9 = 0.243837 11 O19 = 


72.5 = 0.0045021095 v7 = 0.27405120 XK 10-5 
73-6 = —0.059185092 ve = 0.14410406 XK 107% 
-s = —0.0033694729 v = 0.26473430 XK 1075 
76-9 = —0.011665598 vs = 0.12805558 X 10-5 
7x5 =  0.0057683905 v7 = 0.26548789 XK 10% 
75-6 = —0.026336688 vg = 0.12369657 X 1075 
uy = 0.066586701 K 107% ue = 0.057452556 XK 10~% 
us = 0.036619225 XK 1075 us = 0.026845279 XK 1075 


Notice that the average shears are defined as 


1 (3 
Ti} = +) t ds 


The above results for deflections are in close agree- 
ment with accurate values reported in reference 4, in 
which no stresses were reported. 

The same plate is now subjected to a heat loading 
shown in Fig. 4(a) with the stress, o»,)., removed. 
Mathematically, the loading terms, which are of the 
form 

S ad Tds 

are added to the axial edge force-displacement equa- 
tions. For example, to the right-hand sides of the edge 
equations listed above, we add the following: 
edge 4-6: 

h, (1/3) [(aAT)s + 4(aAT)5 + (aAT)¢| 
edge 1-7: 

h, (1/3) [(a@AT), + 4(a@AT)s + (aAT);z] 
edge 1-3: 

h, (1/3) [(a@AT),; + 4(aAT)2 + (aAT)s] 
edge 2-8: 


h, (1/3) [(@AT)2 + 4(aAT)5 + (a@AT)s] 


Similar terms are added to the equations for the other 


edges. Now the solution of the 28 equations is 
o,2 = —18.941310 on = 4.6129141 
o,3 = —88.011599  or5 = —6.2764770 
oy, = 6.7283745  o26 = 8.5620575 
o,¥6 = 25.686173 O27 = —3.2016116 
oys = —7.9721885 or3 = —2.6559770 
oy9 = —14.733092 029 = —34.911030 
Fea-5 = 4.2511075 % = £0.29559233 X 1078 
73-6 = 5.6813984 ve = —0.016503363 X 107% 
75-8 = —3.9719596 % = 0.39932802 X 10-8 
7e-9 = 17.526483 yy = 0.13938995 * 107% 
74-5 = —2.1504006 v7 = #£«0.538908519 XK 10-3 
75-6 = 1.2683979 v = 0.18125297 XK 10-3 
uy = —0.57250006 X 10-* uz, = —0.22470902 X 10-5 
us = —0.30356632 XK 10-3 us = —0.15720397 XK 10-3 


The state of stress in the plate is seen to be complex in 
nature. 
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(II) Cantilever Box-Beam Structure 


Let us consider next the portion of the cantilever 
box-beam structure depicted in Fig. 5. Only one quar- 
ter of the entire structure need be considered because 
of symmetry and antisymmetry. Note that only one 
cover-plate element and only one web-plate element are 
used in the present analysis. The following 12 equa- 
tions are set up to describe the behavior of the struc- 
ture under tip shear load: 


Equilibrium 
cover: 


x = Q: DTo_4 + (1 2)a( oz; + Cn) = 7] 
~Y = 0: a(71-2)¢ — (b/3)(oy2 + 2oys) = 0 


mM, = 0: (b?, 6) (4oy4 — dy2) + 
(a? 6) (2073 a Ors) = () 
web: 
>~Y = 0: (F1-2)0 + (h/4)oy — 76-5 = O 
>~M, = 0: Fa + (th?/24)o,2 — (tah/2)765 = 0 
stringer: 


tai oyA = at| ~ es (To1 Ne + (F214) w } =) 
Force Displacement 
cover: 


shear: (1/2) [724 + (71-2)-]ab/G = 
—(1/2)bu, + (2/3)a(v,; — v3) 


edge 1-2: % = (a/2E)o,» 
edge 1-3: —u = (2/3)boz3 
edge 4-3: v3 = (a/2E) [oy — v(or3 + ors) | 


at wall, zero strain in X direction at point 4: 


Ors — Voy = O 
web: 
shear: (1/4){ [2F/(h/2)t] + (71-2) + 
tes} a(h/2)/G = —(2a/3)x, + (h/2)w, 


Two solutions are determined as follows: for A = 
3.1825, 


oy = 1.8417338 7a-4 = —0.072746390 
o,4 = 1.3082408 Fa». = 0.47294235 
or3 = —0.38977720 7e-5 = £0.79531481 
on = 0.43608027 Fu-2)w = 0.76600623 


u, = 33.079092 % = 368.34675 


v3 = 258.56130 Ww, = 8,551.39836 


and for A = 0 (no flange 


oy2 = 2.9377159 f24 = —0.11603645 

oy = 2.0867510 Fu». = 0.75438169 

o23 = —0.62172648 ts = 0.80112808 

62s =  0.69558366 Fu 2» = 0.75438169 
u, = 52.763854 vv; = 412.42640 


% = 587.54317 w, = 13,141.551 


The deflections, w;, agree closely with the results 
8,740 and 13,350 reported in references 4 and 5, respec- 


tively. In practice, the cover plate should be divided 
into several subplate elements in order to obtain more 
accurate solutions, especially for the stresses. 


(III) Complex Plate Problem 

The plate elements in the previous examples are 
untapered rectangular elements. However, since the 
purpose of the present paper is to present methods for 
general elements, we will consider next the structure 
containing the quadrilateral shown in Fig. 6. The 
structure is assumed to be completely fixed along the 
edge 3-6-9. A single concentrated load of one pound is 
applied at joint 7 in the direction of stringer 7-1. The 
data defining the geometry, gages, and stringer proper- 
ties of the structure are given in Table 1. Notice 
the complex taper pattern of the structure. For pur- 
poses of analysis, the plate element is divided into four 
subplate elements defined by drawing lines connecting 
the bisectors of the sides of the element. It is assumed 
that the resulting lines bisect each other, i.e., Lo; = 
L;-3and Ly, = Ls». This approximation is very accu- 
rate and is considered adequate for practical purposes. 

For purposes of simplicity, the stringers 3-2-1, 
1-4-7, and 7-8-9 are assumed to take axial loading only 
(no bending). Any concentration of load at joint 1 
therefore must be taken by the plate. The following 
69 equations are set up for purposes of solution: 

(1) Equilibrium: three overall equations for the 
plate elements; one each for the six stringer elements 
2-3, 1-2, 1-4, 8-9, 7-8, and 4-7; and Mohr’s circle equa- 


tions as follows: 


4 at joint 4 
4 at joint 6 
2 at joint 7 


3 at joint 3 
3 at joint 9 
4 at joint 5 


4 at joint | 
4 at joint 2 
4 at joint § 


(2) Force displacement: one each for four sub- 
plate elements in shear; one each for six stringer ele- 
ments; one each for four subplate element edges 2-5, 
5-8, 5-6, and 5-4; one each at points 3, 6, and 9 at the 
wall, equating Ke, = o, — vo, = 0; and two at joint 1 
specifying equality of strain in the plate and stringer 
elements 1-2 and 1-4. 

(3) Special: continuity of normal stress at point 6. 

The above equations will not be shown here for lack 
of space. The results are: 


TABLE 1 
Data Used in Complex Plate Problem 


Point x y t Astr Tater. 
1 0.34 1.05 0.06 0.2 0.01 
2 0.17 2.00 0.05 6.7 0.05 
3 0 2.95 0.04 1.0 0.10 
4 0.865 0.525 0.09 0.2 0.01 
5 1.5175 1.7375 0.08 - - 
6 2.17 2.95 0.07 . — 
7 1.39 0 0.10 0.5 0.03 
8 2.865 1.475 0.07 0.7 0.05 
9 4.34 2.95 0.05 1.0 0.10 
3’ 1.085 2.95 0.06 — — 
9’ 3.255 2.95 0.06 — _- 





» = 0.3;E = 10.5 X 108;G = E/2.6 
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On3-6 = —0.27265386 Te-5 = —2.6429064 

On3) = =9= —7.9493600 T-6 = —2.2982464 

One-3 = —4.18371080 m-5 = —6.7777342 

ons =  0.69068216 ra = 0. 51603196 

ons-« = 0.52023626 2 = 0.73893214 

One—-5 = —4.0653728 To-5 = 4. 5846552 

Ons = —2.6686630 m2 = —29.034324 

Ont-5 = 1.97497084 T2-3 = —5.9392944 

Ons-5 = 1.03165332 73-2 = —0.20415276 

Ons-2 = —7.6863910 T35 = 0. 2529434 

One = 3.627194 Te. = 2. 5792232 

Oni-2 = 118.798938 T9-6 = 0.338153E8 

Oni = 119.835898 To-3 = 0. 182082692 

Or} = 106.241466 Ts_9 = 0. 179212648 

Oy i = §51.1507620 T7-3 = —10.319210 

On3-2 = —0.35445004 74= 10319210 

o,9-3 = —0.428864040 ™; = 7.3707550 

Ostr 12 2. 9536096 N14 = 27 . 5453520 

Ostr14 = 1.60556006 n = 41.13978600 

024-1 = —2.2733690 te? = 5.6450818 

Os3s-9 = 1. 2F 864026 3 = 0. 25294344 

o;9-s = 0.676307140 3 = 2.8317380 

Os2-5 = —4.0&60594 ™ = 1. 13668452 

055-2 = —1.72482634 Tr, = —0.48398412 

o,s-5 = 0.22698694 ™é = 2. 5792232 

O;6-5 = —1.37345138 7 = 0 

O,5-6 = —0.67425544 % = 0.62932014 

Os4-5 = —4.2483400 ™ = 0.33815358 

Ty = —2.53115020 

uy = —2.9492236 XK 10-* us = —0.21004946 XK 10-* 
2%), = —0.55663580 X 10% vs = 0.6793866 X 10-6 
U2 = —0.087834044 X 107° uz = —1.64312098 X 10-* 
ve = 0.0509245 X 10* v7 = 1.08154286 < 10-* 
us = —0.62727726 XK 10* ug = —0.24897938 X 10* 


= 1.83670284 X 10~* vg = —0.08168819 X 10% 


U% 

The values of the shear stresses with only one sub- 
script represent values of 7,, at a point. It is seen 
that there is a high concentration of stress in the plate 
at point 1. 

The problem is now repeated taking into account the 
bending of the stringers. The following 30 additional 
equations are set up: 

Equilibrium: two each for six stringer elements, and 
two each for joints 1 and 7. 

Force displacement: two each for six stringer ele- 
ments, and two at joint 7 specifying equality of strain in 
the plate and stringer elements 7-4 and 7-8. 

Also, terms are added to some of the previous equa- 
tions to take into account the effects of the induced 
normal stresses acting on the previously unloaded edges 
of the plate. Again the equations are not listed here. 
The results are: 


on3-6 = 2. 2480124 Os 3-2 2.9224161 
on3’ = —10.538278 Os2-3 = —2.3859745 
One-3 = —1.8791601 Cstri1-2 = —2.7384893 
Cao = 8.3180974 Cstr1-4 = —2.7538833 
on 9-6 = 3.3296684 Os41 = —3.6560925 
Ono—5 = —4.3990600 Oss-9 = —0.14038872 
Ons-6 = —1.4210706 os9-s = 4.3285689 
Ona-5 = —3.5242476 Os2-5 = —3.3834021 
ons—5 = 2.4283620 Oss-2 = 0.51543924 
Ons-2 = —2.7482813 Oss» = —4.3395589 
One 5 = —4.3871115 ose5 = 1.9561518 
On1-2 = 6. 2430413 oss-6 = —0.81177146 
oni = 6. 2548829 Os4-5 = —2.7267520 


5. 1596765 


or 1 


oy1 =  0.21778792 
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= —2.8854254 
= —3.3219767 


—0.35907075 
—0.49090959 
2. 9095313 
—2.7768929 
—2.4879474 

2.3906370 
2.6763490 
—1.0923775 
4.2414952 
—4.8264239 
—5.8253243 
3.3207851 
—3.9891783 


IIb 
Fic. 6. Complex plate problem. 


Ti-~4 
T4-1 
Ti-4 
TI 
T2 
T3 
73" 
74 
TS 
T6 
77 
78 
T9 


‘ 
T9 





3.9891783 


0.080019453 


2.4709443 
3.5661507 
—2.7622167 
—1.0923775 
0. 26365168 
0. 53059266 
2.0842329 
4.2414953 
0.99409506 
0.81520970 
—4.8264239 


= —1.1756672 
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u,; = —0.46135588 x 107% on 3- 
v1 = 0.12099022 x 10-8 és + 
8B, = 0.10642304 x 10-6 on 4- 
ue = —0.17124447 Xx 10-8 On 7- 
ve = 0.013996078 x 10-6 on 7- 
Bo = 0.30975888 x 10-6 On 8- 
u, = —0.5821914 X 1076 Cn oe 
% = 0.30993644 x 10-6 Cote? 

8, = —0.15568449 x 10-6 inne: 
us = —0.13434056 x 10-8 isdn 
vs = —0.12270695 XK 10~® Setr 9 

uz; = —0.64852970 x 107-6 Sete 3-3 
v7 = 0.48160028 x 10-6 Gas 

8; = 0.020781369 x 10~° City 0st 
ug = —0.23011297 K 10-6 

vy = 0.13675498 x 1076 

6s = 0.28212101 x 10-6 

Notice that o, is the stress in the plate element in the 
direction of anedge. The stress in the attached stringer 
iS Os — VOn = Gstr- Note the beneficial effects of in- 
cluding the bending resistance of the stringers. 

Future Work—No mention has been made of the 
bending of plates. This subject requires a separate 
discussion, which will be presented in the future. 

Appendix 
(1) Equation of Shear-Force Displacement 
for an Arbitrary Plate Element 

The equation is derived as follows: Apply a unit 
shear flow, g = 1, around a rectangle that completely 
encloses the plate element as shown in Fig. 1. The 
work done by the induced external loads acting on the 
edges of the plate element is: 

f (6.9 + 6,2) ds (Al) 
where (2) = ( et >») (A2) 
n —sin 20 
Here 6 is the angle measured counterclockwise from the 
x axis to an edge of the plate element. 46, and 6, are 
displacements taken positive in the direction of the 
positive loadings o and 7, as shown in Fig. 1. These 


displacements are expressible in terms of the u and v 
in the form: 


5, ( cos @ sin 60\/u ; 
= ‘ - A3 
(?*) .—sin @ cos )(%) — 
Also, one may write 
ds = dx sec 0 
ds = dy sec 0 — 


After the necessary algebraic operations, the integral 


(Al) reduces to 
fer (u dx — v dy) (A5) 


The internal work done by the unit shear flow is 


simply 


ff (1) Yry dx dy = ff (r:,/G) dx dy (A6) 
Ap Ap 


which leads to Eq. (5). 


MATRIC STRUCTURAL A 


NALYSIS 


11.099617 M, = —0.068812634 
—5.3845391 M. =  0.080470136 
—2.5949072 M; =  0.37518324 
0. 20389746 M, = —0.0037799771 
— 1.3254796 M; = _ 0.14798395 
—1.7708081 Ms, = —0.061923912 
13.981417 My, = —0.027775458 
— 1.3866488 Ri: = 0.17731132 
0.60154132 R,;_; = 0.16750309 
0.3908537 1 R: = 0. 23560983 
0.134138 R; = —0.26041866 
—0.77061277 R, = —0.16826152 
—0.407469 R74 = —0.30077066 
—2.8776203 Ri-3 = —0.30667561 
Rs = —0.10243473 
R, =  0.35583331 


(2) Slope-Equality Condition 

The required equation is derived by passing a pair 
of overlapping cubic parabolas each through four points 
three of which are common to both. The slopes of the 
parabolas are equated at the desired point. Assume 
equally-spaced intervals between points. Let 
parabola pass through the points 7, j, k, and 7’, the other 
through the points 7’, 7’, k, and 7. Equate the slopes 
at pointk. The result is 


one 


fi — 4f, + Of — 4fy + fe = 0 (A7) 
(3) Evaluation of Integrals Involving Axes of Symmetry 
or Antisymmetry 


Let k be located on an axis of symmetry or anti- 
symmetry. Then one may evaluate the integral in 
Eq. (12) as the difference of two integrals as follows: 


Pj Pk Pk 
f f(p) dp = f f(p) dp - { t(p) dp (A8) 
Pi Di Pj 
Then for & on an axis of symmetry, 
Pj 
f f(p) dp = (L/3)(f;r + 4f; + fe) - 
Pi 
(L/3)(f; + 2fe) = (L/3)(fi + 3f; — fe) (AY) 
while for k on an axis of antisymmetry, 
Dj 
f f(p) dp = (L/3)(fi + 4f;) - 
Di 
(L/2)f, = (L/6)(2f; + 5f;) (A10) 


The integrals appearing in Eqs. (12)—(16) have been 
evaluated with the use of second-degree parabolas, as 
have the above results. Presumably, more accurate 
expressions are obtained by use of higher-order parabo- 
las. For example, if k is on an axis of antisymmetry, 
one may derive the following more accurate results: 


Pj 
f f(p) dp = (3L/8)(fi + 2f)) 
Pi 


a (All) 
f f(b) dp = (L/24)(—fi + 14f)) 
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should be used in lieu of the ones derived on the basis 


of second-degree parabolas: 


Pk 
f bf(p) dp = (L?/120) (—f; + 14f; + 47f,) (A12) 
Pj 
Di 
f pf(p) dp = (L7/120)(31f; + 46f; — 17f,) (A13) 
Pj 
Dj 
f bf(p) dp = (L*/240) (23f; + 128f; — 31f,) (A14) 
Pi 
Pj 
f pf(p) dp = (L?/240)(— 3fi + 72f; + 51f,) (A115) 
Pk 
— 


Elastic Instability (Continued from page 858) 
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If k is on an axis of symmetry, the following Eqs. 
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ficant difference is indicated between the critical 
stresses for bending and those for uniform compression. 
Furthermore, in the framework of small deflection 
theory, no significant margin is indicated for the 
critical stress of complicated circumferential stress 
distributions, with respect to uniform compression. 
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Ring With Edge Moments 


Anthony E. Gemma 
Assistant Project Engineer, Connecticut Operactions—CANEL, 
Pratt & Whitney Aircraft, Middletown, Conn. 


May 26, 1960 


gf ox NOTE describes a simple analysis of a ring with edge 
moments: the ring material obeying, in turn, (1) elastic, 
(2) linear viscoelastic, and (3) steady creep (i.e., power creep), 


relationships. 


(1) Etvastic ANALysrts! 
This analysis is essentially from reference 1, and is only in- 
cluded for comparison 
From Fig. 1, the following relationship is found 


r*=r+(r—n)cosg+szsng =~ r+ ze (1) 
and the strain is 


eg = (r* — r)/r = ze/r (2) 


— 


‘sing Hooke’s law for uniaxial stress 
og = Ezg/r (3) 
the radial moment becomes 


h/2 h/2 b 
M, = f oezdzdr = Ee f stds f (dr/r) = 
—h/2 -h/2 a 


¢g( Eh?/12)1n (b/a) (4) 


R = (a + b)/2 (5) 
Then from Fig. 1 
M,Rdé — 2M,(d6/2) 0, M, = MoR (6) 


and the rotation is 
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Fic. 1. 
¢ = (12M6R)/Eh? In (b/a) (7) 
The maximum stresses at r = a at the surfaces z = +h/2 
(o6)mar = +6M)R/h?a In (b/a) (8) 


(2) LINEAR VISCOELASTIC ANALYSIS 


The linear viscoelastic body (Maxwell) for uniaxial stress is 
described by 


é = (1/E)o9 + (1/d)oe (9) 


where dots denote time derivatives. 
From Eqs. (7) and (9) 


09 + (E/A)og = E(z/e)o (10) 


: t ; 
og = e~ (E/M! [ os + Ez rf, ge(E/™! at | (11) 


of = o9 = gE at t=0 


and 


where 


Two cases will be considered: 
(a) If a constant rotation is applied and maintained 
¢=0 (12) 
therefore, 
oo = age (E/Mt (13) 
(b) If a constant rotation rate is applied 
¢ = constant = C (14) 
Then 
09 = e~ (E/Mt { a9? + Aes /r [e(E/Mt - 1} (15) 
(3) NONLINEAR Viscous (OR STEADY CREEP) 
Differentiating Eq. (2) with respect to time 


éo = (2/r)o (16) 
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and using a power law for creep 
ég = Bog", B,n = constants (17) 
we obtain 
og = [(1/B)(z/r)e]!/* (18) 
and the moment becomes 


M, = (¢/B)'/"[2n2/(1 + 2n)(n — 1)\(h/2)A+2)/" x 


(piu — 1) _ qg(*—)) ml (49) 


And the rotation rate is 


¢ = (1/B)(1/MpR)"(h/2)! +2" [(2n2)/(1 + 2n)\(n — 1) X 


(p("- 1)/n 


— gi*-) nm) |" (20) 


REFERENCE 
1 Katz, A. M., Teoria U progosti, pp. 119-121, Moscow, 1956 


——————_+¢ 





Radiation Versus Mass-Transfer Effects 
for Ablating Re-Entry Shields of a Nonlifting 
Satellite 


Ernst Wilhelm Adams 


Aeroballistics Laboratory, Army Ballistic Missile Agency, 
Redstone Arsenal, Ala. 


June 10, 1960 


_ THE VICINITY of the stagnation point of a ballistic satellite 
re-entering the atmosphere of the earth, this note treats 
heat-protection shields consisting of homogeneous, opaque glass. 
Employing evaluations of an exact, transient calculation method, 
derived in reference 1, the note discusses in particular the heat 
balance of the investigated shields (see Table 1) which possess 
different values of the relevant material properties and the sur- 
face radii of curvature, R 

The trajectory data of the vehicle are given in Fig. 1; its re- 
entry angle is 92.35° and the ballistic factor is 0.03 m.*/kg. sec.? 
The material properties of the shields presented in the Table 
include the measured limits of the thermal conductivity k, and 
the surface emissivity ¢ of glasses as functions of the temperature. 
In the expressions employed, u(7) = u* X u,(T) for the viscosity 
and p(T) = A X po, p(T) for the vapor pressure, the functions 
up(T) and p,,p(T) hold for Pyrex glass; u* and A are free, non- 
dimensional level factors. The value A = O eliminates the 
evaporation process, and u* = © yields a hypothetical non- 
melting shield. The constant values used for the specific heat c, 
and the specific weight y~ replace the comparatively small meas- 
ured variation with the temperature of these properties. 

The temperature of the unheated layers of the shields is as- 
sumed to be 300°K. According to explanations in reference 2 
the distance e from the surface to the station 380°K. is used as 
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READERS’ FORUM 





Case R ¥ Cp k X 104 (k/ycp) X 10° p* 
kg kcal kcal m.? 
m m.3 kg. °K m. °K. sec sec 
l 0.635 2,105 0.29 0.61 1.0 1.0 
2 0.635 2,105 0.29 ».00 8.2 1.0 
3 0.635 2,105 0.29 0.61 1.0 1.0 
4 0.635 2,105 0.29 5.00 8.2 1.0 
5) 0.635 2,105 0.29 0.61 1.0 1.0 
6 0.635 2,105 0.29 0.61 1.0 20 
7 230 2,105 0.29 0.61 1.0 1.0 





a scale for the thermal penetration of the shields. The total 


tf 
s=- f v(t) dt 
0 


is given by the integral over the melting speed —v.(1t), extended 
from the beginning of the re-entry at ¢ = 0 to the impact at ¢ = ty. 
According to reference 2, the necessary thickness of the shields 
is'(s + 6) with b = e(t;). The Table presents the following per- 
formance data: b,s,(s + 6), and the maximum surface tempera- 
ture 7;, ,,a,} in addition, the integrals 


1* 
Q= f, q(t) dt 


are presented over (1) the aerodynamic heat-transfer rate gaero(t), 
(2) the amount of heat graa(t) = 1.88 K 107"! K € X T;4 radiated 
back into the air, and (3) the heat convection g(t) within the 
liquid-glass flow. Here, ¢ = ¢* stands for the end of the net heat 
transfer into the shield—i.e., Qaero(¢*) = Qraa(t*). For Case 
1 of the Table, the figure presents some performance data as 


ablation 


functions of the time ¢ 

The Table indicates that the thermal diffusivity k/yc, is the 
dominant material parameter, a result which also holds for the 
IRBM treated in reference 2. A low value of k/ycp yields high 
surface temperatures, steep temperature profiles normal to the 
surface, small melting rates, and hence small values of (s + bd). 
Particularly high surface temperatures, and thus, surface vapori- 
zation rates (obtained from reference 3) and radiation rates, take 
place for k/yc,p = 1077 (see Table). Even then, the blocking of 
the aerodynamic heat transfer Qsero through the evaporation is 
very small as the comparision of Cases 1 and 5 shows. Within 
the total investigated range of the relevant material properties, a 
high percentage of Qaero is radiated back into the air. 

The shields discussed possess a stable performance since, at an 
increase Of Gacro(t) by a factor 8 > 1 (e.g., through a smaller value 
R because of dacro(t) A R*-5), (s + b) rises by a factor 5 with 
1<6<B. 

Even at high surface temperatures, the glass shields investi- 
gated have small values of the thermal penetration, the melting 
rates, and the evaporation rates with correspondingly small en- 
thalpies. Regardless of the temperature limit when disintegra- 
tion of the shield begins at its surface, a thin shield thus may be 
utilized safely up to very high surface temperatures, which con- 
siderably reduce the net heat transfer (Qsero — Qraa) across its 
surface. It is seen that the cooling process of the shields in- 
vestigated is predominantly of a radiative nature with convection 
in the liquid glass and its evaporation playing only minor roles 
in the heat balance. 
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TABLE | 

4 € Ti, maz Qaero Orad d. s 6 s+b 
kcal keal kcal 

"=. m.? m.? m.? mm mm mm 
0 0.8 2,676 69,151 60,194 6,097 5.0 13.8 18.8 
0 0.8 2,489 70,400 43,698 20 , 204 18.1 36.0 54.1 
0 0.4 3,101 67 561 54,258 10,770 7.8 12.8 20.6 
0 0.4 2,753 69,420 32,846 30,559 24.8 34.5 59.3 
1.0 0.8 2,561 61,625 53 , 871 4,805 5.1 13.9 19.0 
0 0.8 2,711 66,245 60,156 0 0 19.4 19.4 
8 2,474 2,359 2.1 15.8 17.9 


49 084 42,795 


Steady Performance Characteristics 
of Gas-Lubricated Journal Bearings 
With Slenderness Ratio L/D = x 


W. A. Gross** 
IBM Research Laboratory, San Jose, Calif 
April 26, 1960 


| Parma AND BURGDORFER! 
characteristics of infinitely 
bearings. These results, together with the approximate side 
flow factors given by Ausman? yield fairly good predictions for 
this note, 


recently published performance 


long gas-lubricated journal 


the behavior of finite gas-lubricated bearings. In 
we give load, friction, and attitude angle results for gas journal- 
bearing films which have slenderness ratio (length/diameter ) 
L/D = wx. These films are long enough to develop pressure 
profiles in the middle of the bearing which are comparable to 
those for infinitely long bearings. They are short enough to per- 
mit a convenient finite-difference solution. The Reynolds equa- 
tion for a steady lubricating film, which is laminar and viscous* 
is 

O(PH%0P/0X)/OX + &PH%0P/dY)/oY AO(PH)/OX (1) 
where X = x/l, Y = y/!, P = p/fa, H = h/h*, with / an appro- 
priate length, and h* an appropriate film thickness. The bear- 
ing number is A = 6yU1/(h**p,). The ambient pressure is pg, 
the viscosity is w, and the tangential surface velocity is U 
Since temperature changes in the film are negligible, the density 
is assumed to be constant in the derivation of Eq. (1) 


** Member of Research Staff, Applied Mechanics 
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BEARING NUMBER, A= 6yUr/(c* pg) 

Fic. 1. Effect of eccentricity ratio and bearing number on load 
W’, friction F’, and attitude angle ¢ for complete journal bearing, 
with L’ = L/D = g-, obtained by digital computer using 
24 X 12 grid points. 
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The film thickness in a journal bearing is given by 
H = ] + ecos B 


where @ is the angular coordinate and h* = c, the radial clearance. 


It is customary to set / y, the journal radius. For computer 
evaluation, the Reynolds differential equation is replaced by a 
difference equation.‘ 

The curves of Fig. 1 were obtained by operating on the pres- 
sure obtained by solving the Reynolds difference equation with 
24 X 12 grids representing perfectly aligned journal-bearing 
films with ambient pressure on the boundaries. Continuity was 
imposed by requiring the pressure and pressure gradient obtained 
from the computer solution to be continuous. The solution to 
the difference equation is accurate to four places, and the accuracy 
of the curves in the figure is about 1 per cent 

The bearing load W is given in dimensionless form by W’ = 

/(DLpa,); the frictional foree F may be obtained from F’ = 
Fce/(uDLU). Note that when F’ 1, we have Petrov friction 
which results when « = 0. The eccentricity and load vectors 
are separated by the attitude angle ¢ 

The attitude-angle curves clearly show the limitations im- 
posed by using a first-order perturbation solution based upon e 
as a parameter. Such perturbation solutions should be reason- 
ably good when e < 0.4. 

It is of particular interest to note the approach of these curves 
to the asymptotic value given by 


(1 + (3/2)e2]!/2f 1 — (1 — 2)!/2 
<a (2) 
€ (1 — ¢*#)*/* 


This theoretical asymptotic value has not previously been veri- 
fied for a finite bearing. 
The friction given is that on the journal 
F’ = Fe/uDLp = J JS [H™ + 3HAdP/d8) dYdB 
ml — aa 24 (3/2)eA 1W’ sin @ (3) 


W’ p—> 


Credit for obtaining the solutions presented in Fig. 1 goes to 
A. Michael, who programed the Reynolds-equation solution, 
and supervised much of the computer effort, R. Togosaki, who 
programed the solutions for load, attitude angle, and friction, and 
I.C. Tang, who ran the programs on an IBM 650 computer and 


collected the results. 
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Bending Stresses Due to Temperature 
in Hollow Circular Plates, Part II 


Marvin Forray and Malcolm Newman 


Design Specialist and Principal Engineer, 
Structures, Applied R & D Division, 
Republic Aviation Corporation, Farmingdale, L.I., N.Y. 


June 10, 1960 


IX A PREVIOUS NOTE,! formulas were presented for the deter- 
mination of the deflection, moment, and shear in a circular 
plate with a concentric hole subject to a linear thermal gradient 
through the thickness. Curves were then given for the cases 
of clamped outer radius coupled with clamped, simple, or free 
inner radius. It is the purpose of this note to present the detailed 
solution with design curves for a hollow circular plate simply 
supported on the outer boundary, and simply supported or 
clamped on the inner boundary. 
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Fic. 1. Nondimensional deflection, moments, and shear for a 
plate with simply supported boundaries. 


(A) Plate With Simple-Support Boundaries 

The formulas for the deflection, bending moments, and shear, 
in nondimensional form, for a simply supported plate with a 
temperature difference Tp between the upper and lower faces, are 


_ wh = 1 J . l1+p a a(i a? r2 , 
aT pb? Ao ) > Ai ~- b b? bt x 
a a? l1+y a r 
k +2 1 — lo l 
al’ | ( “C _ 2 és eas b r 
a? 
| : +: ae log (1) 
b? 
= M, _ 1 , a? 
EaT ph? - 6(1 — v)Ao b? 


- Mo oa 1 ' a? m2 b2 ' Ws 
EaT ph? - 6(1 — v)Ao b? r? "s b | 


Qb 1 a? b? b 
a — + —2 (4) 
EaT ph? 3(1 — v?)Ao |b? a? r 
where 


‘s ste. +4 (22) — 
=a oe ee — 2 oO (oO) 
; 1 + »/\b? a? l1—y rs 


The deflection w is positive downward, M, and Mg are positive 
when the upper fibers are compressed, and a is the coefficient of 
linear expansion. The geometry is shown in Fig. 1 where the 
results are given corresponding to b/a = 0.2, 0.4, 0.6, 0.8, and 
Poisson’s ratio, vy = 0.3. 
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(B) Plate Simply Supported on Outside and Clamped on the 
Inside 

rhe corresponding formulas to Eqs. (1)-(5) for a plate simply 
supported on the outside and clamped on the inside, are 


wh 1 (: _ a a*r? AY ? 
i ~ aee b? he bs og = 
? ia a? i a ,; r ‘ 
2 oe ( - ) (6 
b? b? J b tea b ”) 
M, —1 . ; a? 21 a ' r 
a ~ 26 = 7 ae . oO + 
EaT ph? 12(1 — v)A; \ b? 6 b g a 
e yi a a fh? 41 a? b? (; —v\t 
2 - oO. : - —_ q (7 
\ 7? . b + r? a? b+ pss 0) 
VW —1 2 , 
: — = ——- Jo l1— a 2 log - log : i 
EaT ph? 12(1 — v)A; | 52 b a 
a? _ @ b? a? b? 2a? 
2 +3) log - + eS ee sais 
r? b r? ry? a? b2 


Qo = l [, sa aa b (9) 
EaT ph? 3(1 — v*)A; L b 


where 

wnat (Ho b? 3 + »\ a? . 

~~ 1+ >/ a? 1+) b? 
( bins 4 (1 4)" (10) 
l+y "s b ii ” b —_ 


The quantities defined by Eqs. (6)-(10) are shown in Fig. 2. 
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Fic. 2. Nondimensional deflection, moments, and shear for a 


plate with simply supported outer boundary and clamped inner 
boundary. 
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On Minimum-Weight Rectangular 
Radiating Fins 


Chen-Ya Liu 
Carnegie Institute of Technology, Pittsburgh, Pa 
May 5, 1960 


ARTAS AND SELLERS! have obtained a numerical solution for 
rectangular thin fins where heat is dissipated by radiation 
to surroundings at O°R. In the present note, an exact solution 
for the same problem in terms of well-known functions is pre- 
sented. The fin geometry with least weight can also be deter- 


mined from this solution. 


TEMPERATURE DISTRIBUTION 
The governing equation of the temperature 7(x) in a thin 
rectangular fin is 
d?7T/dx? — (ea/kb)T* = 0 iO < s < ZB) (1) 
where ¢€ is the emissivity of the fin material, o is the Stefan- 
Boltzmann constant, & is the conductivity, and 2) and L are 


the fin width and length. The boundary conditions are 


T = JT atx = 0 (2) 
dT/dx = 0 atx = L (3) 
T=T1 atx =L (4) 


Note that the boundary-value problem is completely defined by 
Eqs. (1), (2), and (3). The additional condition, Eq. (4), is 
merely used as a parameter which will be uniquely determined 

By direct integration of Eq. (1) with boundary conditions, 
Eqs. (3) and (4), the solution for the temperature distribution 
along the fin can be written as 
B(0.3, 0.5) — Birt T)s (0.3, 0.5) = 

(10/3)G ~@/2)(71,/T»)?/21 — x/L) (5) 

where Band B,7,/7,); are the complete and incomplete Beta func- 


tions and 


G = 40kb°/9ecA?T) 


where A = 2bL is the finarea. Eq. (5) gives a solution of one-to- 
one correspondence except when 7, = O which is excluded at 


present and will be discussed at the end of this note. To com- 

plete the solution, it remains to evaluate the parameter 7 ,. 

This can be done by using the remaining condition, Eq. (2). 

From Eq. (5), there results 

B(0.3, 0.5) — Bcrz/T, (0.8, 0.5) = (10/3)G~“/?)( Ty, /To)?/? 
(6) 


After T, is calculated from Eq. (6), the rate of heat conducted 
out of the fin base is given by 


gq = —2kb(dT/dx)x=0 = (1.6eokb)'/2( 75 — T15)'? (7) 


OPTIMIZATION OF FIN 
It is desired to maximize g for a given area A by varying ). 
However, as 6 varies, 7, varies with it according to Eq. (6) 
If Eq. (6) is rewritten as 
p(b, Ti) = B(0.3, 0.5) — B (7r1/7,)(0.3, 0.5) — 
(10/3)G ~"/2)(T,./T,)?/2 = 0 (8) 


where g and are now to be considered as functions of b and 71, 
the problem of optimizing the fin geometry is equivalent to 
finding the extreme values of g(b, T,) subject to the subsidiary 
condition ~(b, T,) = 0. The solution to this problem can be 
obtained by means of Lagrangian multiplier \ defined by 












oo 
J] 
bo 


0g/0b + A(Op/db) = 0 (9) 

0g/O0T, + A(0p/OTL) = O (10) 

Eqs. (8), (9), and (10) serve to determine the stationary values 

of b, Ty, and the constant >». Eliminating \ from Eqs. (9) and 
(10) gives 

T1/T;e = 72~%2(12 — G + G/7(120 + G)!/7}92 (11) 


where the positive sign in front of G!/? is chosen in order to make 
q a true maximum at the stationary values. If Eq. (11) is sub- 
stituted into Eq. (8), there results a transcendental equation 
involving the variable G only. By trial, it is found that for 
optimum fin G = 1.38. 

In the above discussion, the case when 7; = O is excluded. 
For this particular case, the solution of the problem is trivial 
and the optimum condition does not exist because the fin length 
is theoretically infinite. 
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Generalized Integral Relationship 
in Flow Fields With Vortex Sheets+ 


Shinya Kobayakawa* and Lu Ting** 
Polytechnic Institute of Brooklyn, Freeport, N.Y. 
May 10, 1960 


N AN EARLIER PAPER,! a generalized integral relationship re- 

lating the integral of the pressure to the normal component 
of the velocity on the region of a cylindrical surface ahead of a 
Mach plane was derived on the assumption that the disturbance 
potential and its first derivatives in the domain ahead of the 
Mach plane were continuous. It was shown that the integral 
relationship remained valid if, in the domain ahead of the Mach 
plane, there were characteristic surfaces across which the normal 
derivative is discontinuous. ! 

In the present note it is shown that the integral relationship 
is still valid if the region ahead of the Mach plane contains a 
vortex-sheet-like discontinuity surface—i.e., a cylindrical surface 
with generator parallel to the undisturbed flow. Across the vor- 
tex sheet the pressure and normal velocity are continuous. 

The disturbed flow field is divided into upper and lower regions 
by the cylindrical surface composed of the surface of the body 
and wing and the vortex sheet. (There will be more than two 
regions if there is more than one vortex sheet. However, the 
proof can be extended without difficulty. ) 

For the upper region, the integral relationship! is applicable 
and states 


Uv 
f f pnw dS — = f f qn dS + 
Sit Sit 
U 
ff pn@dS — - ff qgrodS =0 (1) 
Set Sot 


wheret S, and Sy are the portion of surface on the wing and body, 
and that on the vortex sheet which lie ahead of the Mach plane 
(X, w) 

x + Biycosw + zsinw) = X 


t+ This research was supported by the USAF through the Air Force Office 
of Scientific Research, ARDC, Contract No. AF 49(638)-445, Project No. 
9781, under the supervision of Professor Antonio Ferri, to whom the authors 
are indebted for his invaluable discussions. 

* Graduate Assistant. 

** Research Associate Professor of Aeronautical Engineering. 

t Notation in this note agrees with that of Ref. 2. 
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and the superscript (+) indicates quantities on the upper region. 
Similarly, for the lower region, we have 


ff n-@ dS “ff dS + 
b n° os = G. 
ames B q ¥ 
Si” Si7 
U 
bn-w dS — = andS = 0 (2) 
t B q 
Sv- 


Sv- 
Across the vortex sheet p and q, are continuous, therefore, the 


values of p n-® and gn on S»* are equal to the corresponding 
values of —p n-W and —gq, on S,~, and the sum of Eas. (1) and 


i 
pnodS = . Qn aS (3) 
B : 
Si 


Si 


(2) yields 


where S; = S,;*+ + S, 
ship for no vortex sheet. Similarly, it can be shown that the 
relationship of the moment of pressure? will not be altered due 
to the presence of the vortex sheet. 

To obtain an approximate solution for the pressure distribution 
on the body surface, which lies in the domain of influence of the 
trailing edge, an alternate form of the integral relationship is 
derived so that the domain of integration is substantially re- 
duced. Let the disturbance pressure p be split into two parts, 
p, and ps, where /p; indicates the disturbance pressure if the vortex 
sheet is replaced by a ‘‘thin sheet,’’ and pf» indicates the correction 
due to the vortex sheet. Determination of p; is equivalent to 
that in the usual interference of wing and body ahead of the 
trailing edge. When the integral relationship is applied to p,* 
and »,~, and then combined with Eq. (3), we have a simple 


Eq. (3) agrees with the integral relation- 


relationship for py»: 


ff pom-@dS = ff (pit — pi~) m-@dS (4) 
Svt 


So* 
where S,* is the surface of the body ahead of the Mach plane 
(X, w) under the influence of the vortex sheet. .S,* is always 


smaller than S;. 
In a manner similar to that employed for the dihedral® and 
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$,*4123+ A137 
Sy=A135 





Fic. 1. Wing and body and the domain of influence of the 
trailing-edge vortex sheet. 
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K\=X-Xo 


Fic. 2. Pressure distribution on the body behind the trailing 
edge 


the wing-body interference problem,’ Eq. (4) can be reduced to a 
line integral which in turn reduces to simultaneous algebraic 
equations corresponding to various sets of values of XY and w 
For the purpose of illustration (Fig. 1), a triangular flat plate 
at an angle of attack @ with a sonic leading edge and an unswept 
trailing edge is mounted as a high-wing configuration on a pris- 
matic body of rectangular cross section, and the pressure dis- 
tribution on the body in the domain of influence of the trailing 


edge is obtained as shown in Fig. 2.5 
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Calculation of Hypersonic Shock-Layer Flow 
Parameters, Knowing Shock-Layer Ratio 
of Specific Heats 


Alfred B. C. Anderson 
Senior Research Specialist, Aerodynamics Div., 
North American Aviation, Inc., Los Angeles, Calif. 


May 20, 1960 


_ RATIO OF SPECIFIC HEATS behind a hypersonic shock, at 
least in the subsonic part, is shown to depend chiefly on 
free-stream normal Mach number Ma and relatively little on 
altitude for a wide variety of engineering-wise flyable trajectories. 
Calculation of hypersonic shock-layer flow parameters is, there- 
fore, facilitated by this ease of picking an appropriate shock- 
layer ratio of specific heats, y2, which is then inserted into simple 
oblique-shock relations, analogous to those of reference 1, de- 
rived for the case where ratios of specific heats on opposite sides 
of shock differ 

Dependence of 7, on Mo (Fig. 2) was determined for two rea- 
sonable extremes of re-entry-vehicle-trajectory design (Fig. 1) 
at Mo = 5, 15, and 25, both behind normal shock at A (Fig. 3) 
and behind the oblique shock at the sonic point. Determinations 
of yo were carried out three ways for equilibrium air; (i) by 


v2 = thr/ho}/{(hr/ho) — (br/Po)(p0/er) (ve — 1)/vel} 


where h = | y/(y — 1)}/(p/p) and where enthalpy ratio h./h,, 


FORUM 3873 


static-pressure ratio p2/p.. and density ratio p2/p.. are appropri- 
ate Mollier-diagram values (references 2 and 3) at the given M 
and vehicle altitude: (ii) by 


v2 = d[In p]/d([In p 


7 


according to Fig. 4 of reference 2; and (iii) by 
£ ‘ 


[2/M,,2 sin? 6] — [ar/a..)] — 1}/}(ar/a.) — 1} 


derived from Eq. (3), below, where the normal-velocity ratio 
tir/i. = pxo/pr is the appropriate Mollier-diagram value 
Circles, dots, and squares are used, respectively, in Fig. 2 to 
represent the three foregoing sets of y2. Curve E, Fig. 2, gives the 
mean of the foregoing ratios of specific heat yz: associated with 
trajectory E, Fig. 1, and curve F, Fig. 2, the mean of the fore- 
going ratios associated with trajectory F, Fig. 1. Results from 
all three ways of calculating y2 agree quite well. Curve M is the 
mean of curves E and F. Fig. 2 shows a very small change of 
yz in going from normal shock on axis to oblique shock at sonic 
point, and, a relatively small change of 7. in going from the high- 
L/D trajectory to the zero-L/D trajectory. For typical calcu- 
lations in and near the subsonic part of the shock layer, one 
might, therefore, pick values of yz associated with curve MV. 

The derivation of the following equations for the different shock 
parameters, for the case where the ratios of specific heats on op- 
posite sides of the shock differ, is based on the equations of refer- 
ence 4 and the following minor approximations for Mo > 1 
and M4 <1: 


[242M — (ve — 1)] * 





27.M,2, [1 + ¥.M@.2]? = y¥.?Mo4, 
[2 + (yo — 1)M,2]/[2 + (v2 — 1)M.2] = (yo — 1 
(y2 — 1)[2 + (v2 — 1)M22] = 2, (1 + yi? 

(1 + yoMe?) = y2/v¥e, [1 + 2yaMa?]/ys = 1 
which holds quite well for air from somewhat below M,, = 5 to 
M,.. > 25. For a normal shock, the Mach number M4, and the 
normal-velocity ratio #i4/%,, are 

Ma? = [2 + (v2 — 1)Mq7)/[272M.? — (v2 — 1)] (1 
and 
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Fic. 1. Dependence of altitude on velocity for high-lift/drag 
and zero-lift/drag trajectories. 

Fic. 2. Dependence of shock-layer ratio of specific heats on 
free-stream Mach number for two foregoing trajectories. Curve 
M is the mean of E and F. 

Fic. 3. Geometry of flow and symbols used 














TABLE 1 





At Sonic Point of Shock 














No OB fia /i7 Pt/Pe aT/Ae PiT/Ptw 
5 66.94 5.24 25.96 2.147 8.86 X 10-2 
67.37 5.24 25.13 2.13 8.74 x 10-2 
15 73.54 11.29 265.2 4.46 7.08 X 10-4 
73.64 11.29 264.8 4.43 7.19 x 10-4 
25 77.22 20.23 787 .0 5.56 6.02 1075 
76.91 20.08 789.0 5.71 6.25 x 1075 

Across Normal Shock on Stagnation Streamline 
Mach 

No fio /tA a4/Q« PA/P PtA/Pto PB/Pea 
5 5. 543 oF 29.68 6.23 X 10-2 17.78 
543 2.23 29.69 6.26 X 1072 17.95 

15 11.8 4.54 289.3 4.59 X 10-4 72.5 
11.8 4.53 289.3 4.58 xX 10-4 166.5 

25 20.58 5.67 830.2 3.81 <x 107% 513.5 
20 58 5.47 833.4 a.7F x ¥O" 508.0 





tia/ihe = po/pa = [2 + (v2 — 1)Mae*]/(v2 + 1)MQ? 
For an oblique shock the foregoing become 


Mr sin? (67 — Br) = [2 + (v2 — 1)M,.? sin? 67] + 
|2y7.M,,2 sin? @r — (y2 —1)] (2) 


and 


r/o = Po/pT = 
[2 + (v2 — 1)M,? sin? 67]/[(v2 + 1)M,.? sin? 67] (8) 

On the same basis, expressions for the velocity-of-sound 
ratio ar/a,, static-pressure ratio pr/p,, stagnation-pressure 
ratio prr/pt~, and Mach number M7 across shock at 7, as well 
as the static-pressure ratio pg/p,. on the body at the point B, 
where the body angle is 6g and the Mach number is Mz, have 
been derived. 

Table 1 compares the values of shock parameters with exact 
Mollier-diagram solutions. Comparisons are made at the point 
on shock where downstream flow is sonic, as well as at point A, 
across normal shock on stagnation streamline. The first line, 
for each Mach number, lists values calculated by the foregoing 
equations; the second by the Mollier-diagram solutions. All 
calculations are based on data in Fig. (1) for the high-L/D tra- 
jectory and the ARDC Standard Atmosphere. 

Agreement in the values obtained by foregoing equations and 
Mollier-diagram-based values is quite close. Similar agreement 
was found for the zero-L/D trajectory, not shown on the Table. 
The difference between the flow parameters for the high- and 


low-L/D trajectories also was found to be small at M,, = 5 and 
15, and relatively small at @, = 25. Most other equations 


in reference 1, not involving Z, have been rederived for case where 
y on opposite sides of shock differs, and have been found to give 
equally good agreement with thermodynamic-chart determina- 
tions. 

The application of the foregoing might be extended to a greater 
range of shock-angle inclinations, more appropriately, by re- 
placing the curves of Fig. (2) with a family, each one of which 
applies to a fixed altitude, and, therefore, a greater range of shock 
inclination. 
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Note on the Sonic Freeze Concept 


J. C. Gibbings 
Department of Fluid Mechanics, University of Liverpool, England 
May 16, 1960 


Ww A STREAM is flowing past a body at just above sonic 
stream speed, there is a detached shock wave upstream 
of it. As the stream speed approaches sonic from above, schlieren 
observation indicates that the shock wave not only moves to 
infinity upstream, but also tends to become plane. Thus the 
upstream boundary conditions at sonic speed are those of uni 
form flow at infinity, whether sonic speed is approached from 
above or below; they have no discontinuous change. 

As the Mach number upstream of the shock, 14, approaches 
1.0 from above, then, because the shock is tending to a plane one, 
so also will the Mach number downstream of the shock, J/., 
approach 1.0 and 


(M2? — 1)/(M\? — 1) > -1 
or dM, = —d., (1) 


Now the Mach number MM at a point on the body, for subsonic 
free-stream flow of Mach number J), can be written 


M = f(M ) (2) 


ub 


and the temptation here is to identify MW .u». with M2 for small 
but finite values of 61/, (see references 1, 2, and 3). This, how 
ever, leads to the variation sketched in Fig. 1, whereas the ob 
served experimental variation is as sketched in Fig. 2 

This difficulty can be avoided. From Eq. (2) 


4M = fow.'dM, 


where f’ denotes df/dM,. Then at sonic stream conditions 


dMwv.* = feav.’* dM (3) 


where dM.y.* and feup.’* are limiting values at M) — 1.0 from 
below. 

From the previous discussion of the continuity of the upstream 
boundary conditions at sonic speed, it is not unreasonable to 
assume that there is no discontinuity in f’ at sonic speed. Thus, 
on approaching sonic speed from above 


AMyuyp.* = fav.’* dM2 


= —fay.’* dMo 


from Eq. (1). Comparing this with Eq. (3), one has 


(dM/dMo)suv.* = —(dM/dMo)sup.* 


M 
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and as a discontinuity in f’ has already been excluded 
(dM/dM,)* = 0 


This relation, first obtained by Liepmann and Bryson,’ is 
satisfied by the experimental results typified in Fig. 2. It is 
here derived without the necessity of identifying M, with MJ, 
which leads to the incorrect variation of Fig. 1 
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On Almost-Free-Molecule Flow Through an 
Orifice7 


V.C. Liu and G. R. Inger** 


Professor and Graduate Student, Respectively, 
Department of Aeronautical and Astronautical Engineering, 
University of Michigan, Ann Arbor, Mich. 


May 20, 1960 


{ bew APPLICATION of the principle of almost-free-molecule 
low, which is essentially a form of first-order Knudsen 
iteration of rarefied-gas dynamics, has shown some very en- 
The effusion of rarefied gases through an 


2 


couraging results.!: ? 
orifice into a vacuum is a very instructive problem for the pur 
pose of studying the basic nature of rarefied-flow phenomena. 
The object of the present analysis is to provide a microscopic 
description of the flow parameters pertaining to the steady 
effusion from an orifice, the diameter (D) of which is of the same 
order or smaller than the mean free path (A) of the reservoir gas. 
A thin diaphragm (t/D << 1) which has a small circular orifice 
separates a large high-pressure (p,) reservoir from the low-pres- 
sure (p2) region. The pressure ratio will be assumed large enough 
(pi/p2 > 10%) to permit neglect of the back flow from the low- 
pressure side. This condition distinguishes the present problem 
from the pitot-pressure problem of reference 1. 

Consider first the case of free-molecule effusion, where \ >> D 
and molecules move through the orifice essentially independent 
of each other. The deviations of the resulting molecular dis- 
tribution from its equilibrium state will be negligibly small and 
promptly wiped out by the intermolecular collisions, which always 
tend to set up and preserve the equilibrium state. The loss of 
molecules through the orifice, however, develops a trace of mass 
motion toward the orifice due to absence of those collisions that 
the lost molecules would have made with the ambient molecules 
This trace of orifice-bound mass 
The principle 


on their return from the wall. 
motion grows in prominence as A/D decreases. 
of almost-free-molecule flow is applied here to calculate the 
molecular flux of the mass motion within the reservoir as a result 
of intermolecular collisions or their absence. The following 
physical model is devised for this purpose. 

Imagine the orifice were closed with an imaginary disc of diam- 
eter D as the orifice; then equilibrium (Maxwellian) distribution 
of molecules in the reservoir would be restored through scattering 
of reflected molecules, from the imaginary disc, with the ambient 
molecules. It is postulated that the net molecular flux toward 
the orifice, inhibited by the scattering action of the imaginary 
reflected molecules, is equal to the difference between the true 


t This investigation was part of a broad upper-air research program 
supported by the USAF Cambridge Research Center under contract No 
AF19(604)-5477 with the University of Michigan 

** Now with Aerodynamic Research Group, Missile and Space Division. 
Douglas Aircraft Co., Santa Monica Calif. 
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Fic. 1. Geometry of scattering analysis. 
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effusion flux through the orifice and its calculated value based 
on free-molecule hypothesis. It is assumed that this molecular 
flux, due to absence of scattering, amounts to only a small frac- 
tion of the corresponding free-molecule effusion flux. This im- 
plies that the present theory is valid only when D is not much 
larger than XA, so that the single-collision analysis used here is 
acceptable. The rate of collisions (N ja) between the molecular 
rays incident on, and reflected from the imaginary disc is calcu- 
lated for rigid spheres on the basis of classical kinetic theory 
The cross section of the spheres is taken from the experimental 
determination of the mean free path. It is assumed that every 
collision throws the incident molecule out of the impinging 
stream on the disc. This assumption can be justified only when 
\ is not very small compared to D. 
“‘closed”’ reservoir are now in thermodynamic equi- 
To determine the net 


We assume that the molecules 
inside the 
librium with Maxwellian 
contribution, through intermolecular collisions, of the hypo- 
thetical reflected molecules, we need to calculate the molecular 
flux (Naa), originated from the diaphragm surrounding the 
orifice, that is thrown into the incident ray through their collisions 
These 


distribution. 


with the reflected molecules from the disc molecules 
(Naa) would otherwise not impinge on the orifice area. 


the net impingement inhibition for the imaginary reflection ray 


Thus, 


is equal to (Nia — Naa). 

Let Na be reflected molecular flux from the disc, which is equal 
to rD*nc/16 where nm = molecular density, c = mean thermal 
velocity; E(E = 0.92/d for air) the average collision expectancy 
per unit distance traveled by a single molecule moving through 


molecules under equilibrium (Maxwellian) distribution.* 
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Referring to the coordinate system shown in Fig. 1, we obtain 


SCIENCES NOVEMBER, 1960 


Nia E 2r 2r , Cg ee D/2 , x polo sin 0) cos? Oe 2ry/d , 
a dd dy A dpy : —- a | 
Na tr? 0 “ 0 : 0 0 . 0 [ro + po? + 270 po sin O cos ( yo — 0)}3'? 


Let p and ¢ be the polar coordinates of an area element of the diaphragm outside the orifice area; r the distance from the point rp to this 


elemental area; @, the angle between 7 and the normal to this elemental area (see Fig. 1). 


We obtain 


Nad E - - 7 wits ‘ a © pporo? sin By cos? Aye~ 7%o/d 
a dd» dy do d% dp dpy aren di 2 
ie or tO 0 0 0 D/2 0 0 ve 


After the exponential functions in the integrands are approxi- 
mated with appropriate parabolic representations (these approxi- 
mations remain close when D is not much larger than A), we ob- 
tain from Eqs. (1) and (2) 


Nia/Na ™& 0.153 (D/d) Fy(X/D) (3) 
Naa/Na & 0.019 (D/X) Fo(X/D) (4) 


where F;(A/D) and F.(A/D) are given in Fig. 2. 

From the hypothesis that the intermolecular-collision effect 
on the effusion rate is equal to the molecular flux inhibited from 
hitting the orifice area when closed—namely (Nig — Naa), we 
have 


N/Ner =1+ [(Nia — Naa)/Nr] (5) 


where Ny is the free-molecule effusion rate through the orifice 
(note that Ng = N,» as a first approximation). A graph of 
N/Nr as a function of \/D is shown in Fig. 3, in which are also 
shown the experimental results* which were taken with p;/p2 > 
10%. Itis felt that with refined analysis—e.g., taking into account 
the contribution to the impinging flux from molecules, emerging 
from collisions (N;), that are deflected toward the orifice area, 
the validity of the theory can be extended to lower values of \/D 
The present analysis appears fruitful in illustrating the functional 
dependency of the effusion-flow rate on the intermolecular 


collisions. 
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An Equivalence Principle 
for Water-Exit and -Entry Problems} 


John P. Moran* 

Graduate School of Aeronautical Engineering, 
Cornell University 

June 3, 1960 


5 beers EXISTS an extensive body of literature concerned with 
the impact forces exerted on a body entering water, a review 
of which was recently made by Szebehely,! while a smaller, but 
growing, number of analyses deals with problems of water exit.2* 
However, it does not seem to have been recognized in these pre- 
vious studies that, outside of viscous effects, there is a definite 
equivalence between problems of water-exit and -entry (provided 
that there is no cavitation). 


+ Based on a portion of reference 8 in which the principle was demon 
strated for slender bodies. The author is indebted to Dr. W. R. Sears for 
suggesting the problem and supervising the research, and also to Dr. S 
H. Lam, for pointing out that the existence of the principle does not depend 
on the slenderness of the body. 

* Formerly, Graduate Student. Now Associate Research Scientist at 
THERM, Inc., Advanced Research Division, Ithaca, N.Y 
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Fic. 3. Effusion rate vs. \/D 





To simplify the demonstration of this equivalence principle, 
we shall consider the vertical exit or entry of a symmetric body 
in uniform axial motion, as shown in Fig. 1 The concept is 
readily generalized to other situations, as will be noted below. 

As usual, we assume the flow to be incompressible and irrota 
tional. We may, therefore, work in terms of a potential ¢g, whose 
gradient is the fluid velocity and which satisfies Laplace’s equa- 
tion by continuity: 


2 = 0 (1 


An integration of the Eulerian equations of motion yields Ber- 
noulli’s equation, which in the space-fixed coordinates (x’, y’, t’) 
defined in Fig. 1 is 
Oe ot’ + (1/2)(Ve)? + (p — ps)/p + gy’ = 0 (2) 

where p denotes the pressure, g the acceleration due to gravity, 
and the constant of integration has been evaluated on the water 
surface, y’ = y,’, far from the body creating the motion 

Neglecting the density of the air above the surface relative 
to that of the water below, we require that the pressure on the 
surface, ps, be a constant. Thus, from Eq. (2) 


gyvs’ = — [d0¢/dt’ + (1/2)(Ve)? ] y= ye’ (3) 
It is somewhat more convenient to work in the dimensionless 


body-fixed coordinate system (x, y, Z) defined in Fig. 1 and re- 
lated to (x’, y’, t’) by 


x=x’/L, y= Z—y'/L, Z= +Ut/L, 


d(x, y,Z) = g(x’, y’, t')/UL (4) 

Where a double sign is indicated, the upper sign refers to the 

exit problem and the lower to the entry situation. The free- 
surface boundary condition, Eq. (3), now transforms to 

ry,//U2 = —[- is 9 2 (5 

gys /l [+Do/Dt + (1/2)(Vo¢) er ) 


where D(_ )/Dtdenoteso( )/dy +  )/dZ. 
The boundary condition of no flow through the body may be 
written 
P , 06/O0x 
dX(y)/dy = + : (6 
U + 06/dy _Jx = X (y) 


Now taking note of the homogeneity of Laplace’s equation, 
Eq. (1), and the way in which the signs change in the boundary 
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conditions, Eqs. (5) and (6), according to the direction of motion, 
we see that the solutions for the water-exit and -entry problems 


can be unified if 
@ = +(x, y, Z) (7) 


That is, the potential of the flow about a body leaving the water 
differs only by a change in sign from that of the flow about the 
same body in entry 

Now the pressure formula, Eq. (2), may be made dimension- 
less and put into the body-fixed coordinates by substitution of 
Eqs. (4) and (7) 


C(x, y, Z) = 2(y — Z)Fr7! — (V)? — 2D@/Dt_ (8) 


I 

where Fr is the Froude number, U?/gL. It is thus seen that the 
pressures felt on a given body, passing through the water surface 
at a given speed, depend only on its location and orientation with 
respect to the surface, but not on the direction of its motion. Fur- 
ther, since in the inviscid approximation, the forces felt by a body 
are obtained by integration of the pressures, it is seen that both 
the magnitude and direction of the axial thrust are, in this case of 
vertical surface crossing, similarly independent of the direction 
of motion. Since the viscous drag always opposes the direction 
of motion, however, this equivalence principle cannot hold in a 
real fluid. But if viscous effects are either small or predictable 
theoretically, the principle permits a simplification of the rather 
unwieldy experimental procedures required to study water-exit 
problems.? One could instead perform the easier entry experi- 
ment, relating the results afterward to an equivalent exit situ- 
ation 

As noted earlier, it has been assumed in this analysis that 
cavitation is not present. Should a cavity form, it would 
destroy the symmetry of the boundary-value problem posed by 
Eqs. (1), (5), and (6). The equivalence principle is thus limited 
in its practical application to the surface crossing of slender 
bodies at relatively low speeds. For example, theory predicts 
that a biconvex airfoil of 30 ft. chord in vertical passage will 
cavitate at about 150 [t./sec. if its thickness ratio 7 is 5 per cent, 
and at 100 ft./sec. if r = 10 per cent,’ assuming that the critical 
pressure for the inception of cavitation is the vapor pressure of 
the fluid 

Note that the derivation performed above is equally applicable 
to plane and axisymmetric bodies; in the former case, X(y) is 
the ordinate of a symmetric airfoil, while in the latter, it is the 
radius of a body of revolution. Finally, note that the equi- 
valence principle applies only to the vertical component of the 
forces felt by a body, insofar as the sense of the force is concerned. 
The direction of the horizontal component depends on the direc- 
tion of the horizontal motion, as is immediately evident if we 
consider a motion parallel to the free surface. 


Ay 


*U 






FREE 
SURFACE 
X(y) 


'? 
Fic. 1. Symmetric body in vertical water-exit. 0,0’ origins 
of body-fixed and space-fixed coordinates, respectively. X, x, 
y, Z dimensionless. normalized by L. 
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Dynamic Response of Beams With Large 
Amplitudes 


Thein Weh 
Southwest Research Institute, San Antonio, Texas. 
June 6, 1960 


HE DYNAMIC RESPONSE OF BEAMS with large amplitudes has 
been studied before by several authors.! However, these 
studies are somewhat limited in that they assume simply sup- 
ported ends or are confined to a study of either free vibrations or 
the steady-state motion, the transients being ignored. This is 
not surprising, since, owing to the nonlinearity of the equations, 
superposition of free and forced vibrations is not possible 

The method outlined below is free of most of these restric- 
tions. It is emphasized that the method can, in principle, be ap- 
plied to clamped-beams although the calculations are consider- 
ably more tedious than in the case of simple supports. For 
the purposes of this brief note, simple supports are assumed, 
and the forcing function is taken as an exponential 

The equations of motion of a beam of unit width with immov- 


able edges and large deflections may be written 


Eh f. 1 =) 
= — dt 
2 J0 \o 


(1) 
Ow 12Na? o*7 O77 
Se ae ee “Tie 2 = goe ** 
ot! Eh® os? or? 
where 
t=x/a, @ =w/a, ¢ = (ht/a*?) V(E/12p 


N is the axial tension, E is the elastic tensile modulus, a the length, 
h the depth of the beam, w its lateral deflection, x the linear 
coordinate, p the mass density of the beam material, and ¢ is the 
time. 

If the physically plausible assumption is made that the tension 
N remains constant at its initial value during a sufficiently short 
interval of time A, then it is readily seen that the second of Eqs. 
(1) assumes a linear form and may be integrated. By consider- 
ing successive intervals A, the nonlinear problem may be solved. 

Consider a simply supported beam, and assume that during a 
short interval A, N remains constant at N. 


Let & = =T,(¢) sin mw= (mn Odd) (2) 


Expanding the forcing function in terms of the normal modes 
and substituting Eq. (2) into the second of Eqs. (1), one obtains 
a differential equation for the time function 7,. The solution 
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of this equation may be written 





A ‘ ie of the initial displacement and velocities, respectively, in the nth 
T, = - = Gn} cos prt + + Qn} sin prt + qne mode, and p, is a nondimensional frequency given by 
a apn Pn 
(3) n? = n?x?2[n?x? + N(12a?/Eh*)] 4 
Remembering that these equations are valid only for the interval A, and that N is the value of the tension at ¢ = 0, substitution of Eqs 


(2) and (3) into the first of Eqs. (1) yields 


2, n'*x*Ek A,?* 


N= (n = odd) (5 
2 4 a? 
n=1 
Elimination of NV between Eqs. (4) and (5) results in 
pr? = ne’) + (3/n*h?) zz mids (n, m odd) (6 
m 


If the final displacements and velocities at the end of the interval A are made the initial conditions for the second interval A, and so on, 
and time is reckoned anew at the beginning of each interval, one may write the following equations for the quantities A, v and p at the 
beginning of the (ry + 1)* interval 


An, 7. 1 i Uny 1 a ° —ad 
= — Qn, } COS (Pn, ,Ar-1) + + da; sin (pn,, 4-1) + Gn, 6 "=! 
a a apn, Pnr-1 
(j 
Uny A nr—1 . Vnr—1 —aA 
= — Pn, — Gnp—,}) sim (Pnr_, Arp—1) + + Gig, @ 9 COB (Py, Ay.) — Cay € 
a a a 


and 





3 = 
— res 24 2 2 
Pn,-? = n‘x4 (: a nth? > mA m, ) (n, m odd) (8) 


m 

Eqs. (7) and (8) can be applied repeatedly to the problem for any 3 = f, fa — f') dn = 0.65 
odd integral value of m and the response of the beam determined 
at any time by Eqs. (2) and (3), until the desired accuracy is 
attained. It is particularly noted that this procedure yields not 
merely the steady-state solution but the transient response as 
well. 

It is necessary to add, however, that references 1, 2, and 3 are, 
within the limitations noted in the first paragraph, more ‘‘exact’’ 
mathematically than is the present method, which is an approxi- 


f(n) being the usual Blasius function. This technique was pre- 
viously used by Rott and Hartunian? to estimate the recovery 
factor corresponding to their interpolation for velocity profiles 
in the boundary layer behind a moving shock; unfortunately, its 
application to the exact profiles for this flow, although straight- 
forward, does not lead to especially simple results, on account of 
the dependence of the function f() on a shock-speed parameter, 
so this application is not discussed here. 

If we treat pu, the product of density and viscosity, as constant, 
then the Blasius solution for the flat-plate boundary layer holds 


mate numerical procedure. 





REFERENCES : ible fl If v is th f ; 
: ee . -ssible yr. s the stre: : so the 
Woinowsky-Krieger, S., The Effect of an Axial Force on the Vibration of eee eee eee ¥ is the stream function, so that 
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2 Burgreen, David, Free Vibrations of a Pin-Ended Column With Constant 
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June, 1951. ia 
3 Eringen, A. C., On the Nonlinear Vibration of Elastic Bars, Quart. Appl. — 1/2 f/ <x Ua y p 
Math., Vol. 10, p. 361, 1952. = emer" Fah ¢@ te o \p. dy (1) 
where suffix ~ refers to the constant conditions in the free stream, 
we find 
+ u/us = f'(n), 1/tw = f"(n)/f"(0) (2) 


where 7 is the shear stress u(Ou/Oy), and 7, its value at the wall. 
The equation of motion becomes the Blasius equation, which 
may be written 


A Note on the Recovery and Reynolds-Analogy 
Factors in Laminar Flat-Plate Flow 


f= —f'"/f" = —(d/dn)(log f”) (3) 
D. A. Spence* with boundary conditions f = f’ = Oat 7 = 0, f’ ~lasy7— o. 
Cornell University, Ithaca, N.Y. The distribution of enthalpy h = c,T across the layer is given by 
June 9, 1960 Crocco’s integral of the energy equation: 

¥ dh x {r\o-! 
[ IS WELL KNOWN from the computations of Pohlhausen, h-hh, = du — 
Crocco, Emmons and others, described, for instance by du} w Te 

Young,' that the recovery factor for the compressible laminar edt ¢ yen sfe\i-e 
boundary layer on a flat plate is approximately ¢!/? for Prandtl 7 Jo ea du e i, du (4) 


numbers o near unity, and the Reynolds-analogy factor 2c;/cy is 
approximately o~?/3. The purpose of this note is to point out a 
simple way of deriving these results by expanding the Taylor 


From this expression, after substitution from Eq. (2), the re- 
covery and Reynolds-analogy factors are found to be 


series for ry and s in powers of — 1); in this way it will be , , , es eo 
eries torrvra : p (@ 3 t y it ill ai (hy = ho)/[(1/2)ux?] = 2¢ (f")? dy (f")2-¢ dy 
shown, analytically, that r/o!/? is 1 + O(@ — 1)*. The same 0 > 
form is found for so”, where (5 
egg: ~ r @ '’ g oni 

* On leave from the Royal Aircraft Establishment, Farnborough, Hants, s = 2 ‘cy = ca fi a dn (6) 
England. 0 a 


in which gn = [490/(a? + pn?)nw] and A, and v, are the amplitudes 
















in 


sé 


(G 





nplitudes 
n the nth 


1 of Eqs 


(6) 


d so on, 


b at the 


- 
| 


iS pre- 
over) 
rofiles 
ly, its 
aight- 
unt of 
neter, 


stant, 
holds 


(1) 
eam, 


(2) 


wall. 
hich 





READERS’ 


where a = f"(0). Here h, is the wall enthalpy at zero heat trans- 
fer, c, is the Stanton number 


(R/Cp)(Oh/OY)w/ pothal hy — hw) 


and 


Ce = Tw/(1/2)pctta* 


Clearly y = s = 1 when o = 1, and we now differentiate the ex- 
pressions with regard to o to find their variations near this value. 


Recovery factor: 
From Eq. (5) 


d i r 9 vf f”)? d: 
tn a? 


) "gry 2— ” "\y fr\2— { we 
Joss) f, cry van - f, (log f”)( f”)? * dn¢ (7) 


Integrating by parts and using Eq. (3), the right-hand side be- 


— Qgi/2 f, (ry dn Pg dn r Cid o dn (8) 
d r 1 2 fr dn f," sya 
de he's ar ne Poa oe 0” ” 0 d ] 


comes 


and 


= f dad _— ping f’ ” 1 
E | ‘ f" dy 
1 | y 1 
=_— - f’2 = = - 
2 0 2 


ie. (d/do)(r/o'?)¢-1 = 0. The second derivative may be found 
in the same way. Differentiating Eq. (8) 


d 1\2/ o — _ ee , 
clogs") fy'san fy f")?—¢ dn - [han foes cs)? anh 


Setting o = 1 
alg fh 3/4 2f, f" dn 2(1/2)f%(log f” 
da? \aV?},=1 si ae 2 (1/2)f(log f°) 
| ig; ‘(log f”) f, sr d 
o J] Slog f") + Jy Sf dn Ian 
. J £ ( 
— ” | f3 : 
2. sary f. Pang 
2 f, fl —f’) dn (9) 
This integral has been evaluated as 0.73096 from a ‘‘marching’’ 
solution of Eq. (3), starting from the accepted value a = 0.33206 


V2, at the Cornell Computing Center. With this value 
(d?/da?)(r/o"?),—, = 0.01904, and as stated 


r/o’? = 1+ 0.0095(o — 1)? +... (10) 


Reynolds-anology factor: 
From Eq. (6) 


d 1y\ 1 od ‘oe Filial 
> = = = - 11 
(2 ‘) Ss - f. (108 a \E) ‘i ats 
and 
d 1 co 1 ss mn oF : 
Px (:) -—- l= i (tos ys dyn = — f, fa —f')dn 


(12) 


The same method of evaluation gives 
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f, fl — f’) dy = 0.35115 


whence 

(d/do)(1/s)em1 = 0.64885 = n (say) 
In this case, an alternative method of evaluation is possible from 
the asymptotic expansion used by Blasius himself (see, for ex- 
ample, Schlichting’, p. 105) since 


fi fl —f dy = - 
2)f?} 


[log a + (1 2)f?] = lim Slog (a@ ‘f”) — (1/2 j 


0 eo 
and for large n 
f” ~ vy exp [—(1/2)f?] 

The required limit is therefore log (a/v), where, because of the 
factor ~/2 used in defining », a and y are 1/2 times the values 
0.332 and 0.231 given by Schlichting; thus, log (a/y) = 0.3628 
This differs by nearly 3 per cent, a surprisingly large amount, from 
the computed value, but the latter is more likely to be accurate. 

The second derivative of 1/s with respect to o is found by 
differentiating Eq. (11): 


d_ 1\81 _ ay 
da oJ) s a 
f”\2 f" go 
(oe: y¢ ) dn 
Qa Qa 
} 3 (1oe)' dn 
—2 4 1—f’)x| 
So ial (13) 
"iad 
f (ioe) dn 
a 
2f, f dn X 
f fl —f') dy 
” | 


= A (say) } 


, a 
| % 
i) 

| 

bo 
S\ a 
eee” 
— 
maim 
a 
iT 

+ 

ll 


ll 


integrating by parts. The value of A has been computed* as 
0.48879, and writing n = (d/da)(1/s)g.1 = 0.64885 
(d?/da?)(1/so")gx) = A — n? + 3n — 1 = 0.69421 


Therefore, as stated 
s=oa"[1 — 0.347(¢ — 1)? +...] 
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* Alternatively, A can be found from a previous result, since it equals 
So” f1 — f")dn + [log (a/y)]? = 0.4972 


but again, the computed value is probably to be preferred. 





Effects of Flow Unsteadiness in Turbomachines 


G. T. Csanady 
Senior Lecturer in Mechanical Engineering, 
University of New South Wales, Australia 


June 14, 1960 


\y IMPROVEMENT on the classical mean-streamline theory of 
energy exchange in turbomachines is to assume statistically 
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steady flow through the rotor. Consider an axisymmetric 
control surface wrapped around the rotor; the total external 
torque M (shaft torque plus friction torque on the control sur- 
face), when averaged over a long enough period to even out 


fluctuations, becomes 
M = Sf (pC.C. 1) dS (1) 


where the symbol (_) means a time-average, the integral extends 


over the whole control surface .S, and 
p = density 
C,; = tangential velocity component 
C, = velocity component normal to surface element dS, to 
be counted positive if pointing outwards 
y = radius from axis of rotation 


In the incompressible case p = const. but the velocity com- 
ponents fluctuate and may be resolved into time-mean and 


fluctuation 
C; = (C,) + C,’ 
=(C) 46, 
so that the mean-velocity product becomes 
(C.C,) = (CiMKC.) + Cr’C, (3) 
i.e., it involves a velocity correlation. Hence, the torque may 
be written from Eq. (1) as 
M= Sf AC) (C) rdS + Sf ACI'C,’) r dS (4) 
The quantity p(C,’C,’) has the physical meaning of a shear 
stress, analogous to the Reynolds stresses of turbulence, and the 
second term on the right of Eq. (4) gives the torque caused by 
these shear stresses. The first term is the conventional mean- 


momentum-flow torque. 

When the fluid is compressible, the density also fluctuates and 
one obtains the further correlation-torque contributions 

S (0'Cr’XC,)r 1S + S(p'C.’ Cdr dS + f(p'Cu'C.')r dS 

In the presence of density gradients, terms like the above will 
cause mass transport. 

From an energy balance, one may also write down an ex- 
pression for power (shaft power less friction on the control sur- 
face). The energy transport outward is 


J u+ i pC, dS 


where u = internal energy per unit mass, C = velocity magni- 
tude. The work done by pressure forces on the boundaries per 


S pc, dS 


Hence the mean power input: 


/ C2 
P= f ( (ou + p = +p}C. )dS 


In the incompressible case, the above reduces to 


unit time is 


\ 


/ C2 \ 
P= Sl 2 + pic, jas 
i \ 
/ C2 
= f (>( : \ + (p) 
hi, es 
of ( > )«C.) dS + fiver as + 
\ & / 
eae J . rio ih 
p (C'C.’) (C) dS + of - Yds 
Fa y 


(6) 
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Of the correlation terms on the right, the first represents trans- 
port of ‘‘turbulent’’ energy, the second, a pumping action of pres- 
sure fluctuations on the boundary. When the density also 
fluctuates, a number of further terms appear. 

The fact that velocity fluctuations may sustain a shear stress 
in the fluid offers a possible explanation for the drop of pump 
head with distance from impeller periphery observed by Acosta 
and Bowermann.! It may also have a bearing on the problem of 


prerotation. 


Examples 


(a) An isolated, straight cascade in motion in a perfect fluid 
Let velocity components normal to cascade front (x — axis) be 
u, parallel to cascade front v. The velocity correlation of 
interest is (u’v’) and the time-average can clearly be replaced 
by a space-average along a line parallel to cascade front, over 
one blade spacing ¢, in a frame of reference fixed to the blades 


(4) 


The velocity-fluctuations in this expression are deviations from 


the ‘‘space-mean”’ velocities 


l t 
udy, (v 
t J0 , 


By continuity, the mean axial-velocity component (4) equals 
the axial velocity at upstream and downstream infinity, respec- 
tively, ahead of and behind the cascade. Also, by Kelvin’s 
theorem applied to control circuits upstream and downstream, 
extending to infinity, it follows that the mean tangential veloci- 
ties are equal to their values at infinity. Considering, then, the 
equilibrium of the fluid within control circuits upstream and 
downstream, respectively, it follows that the correlation stress, 
as given by Eq. (7), must be constant with distance from the 
cascade, because the force due to mean-momentum flow is zero 
The value of the correlation stress is zero at infinity, hence it is 
zero everywhere. This conclusion may also be obtained from a 
solution of the equation for the velocity potential upstream of a 
cascade. ? 

(b) Two cascades in relative motion. The problem cannot 
now be reduced to a steady-flow one by fixing coordinates on one 
of the cascades. However, by Kelvin’s theorem applied to a 
control circuit, which includes a large number of blade-spacings, 
it still follows that the correlation stress is constant with dis- 
tance from the cascade. This result applies separately: (1) in 
front of the first cascade; (2) between cascades; and (3) behind 
the second cascade. Clearly, in regions (1) and (3), the correla- 
tion stress will again be zero, but it may be finite between the 
cascades, giving rise to a kind of ‘‘blade interference.’’ Apart 
from two successive blade rows, similar interference may exist- 
e.g., between the impeller and volute of a centrifugal machine. 

In a viscous fluid one finds that the swum of viscous stress and 
correlation stress is constant (e.g., zero ahead of the cascade) 


with distance from the cascade. 
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